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The  spatial  and  frequency  distribution  of  path-corrected  Rayleigh 
waves  from  the  San  Fernando  earthquake  arc  systematically  related  to  the 
faulting  process  no  n  multiple  event.  The  surface  wave  source  in  taken 
to  be  a  depth-distributed  set  of  double  couplca.  A  least-squares  in¬ 
version  is  used  to  find  the  net  of  source  parameters  which  optimally  fit 
the  variance-weighted  data.  In  thic  study  the  strength,  the  clip  angle, 
and  th..  rupture  delay  times  arc  optimally  determined  for  the  spatially 
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20.  Abstract  Cont’d. 


assumed  distribution  of  double  couplea.  A  sophisticated  error  analyst! 
is  performed  to  estimate  the  uncertainties  oi  the  calculated  model 
variables. 

One  of  the  interesting  side  results  from  this  inversion  technique 
is  the  a  priori  determination  of  what  azimuthal  distribution  of  obser¬ 
vations  is  nost  beneficial  in  determining  the  source  parameters,  the 
uncertainties  in  the  calculated  model  parameters,  and  the  relative 
coupling  between  various  parameters.  Of  the  18  ITh’SSN  stations  used  for 
this  event,  two  stations  contain  a  total  of  30 7.  of  the  total  information 
in  the  data  set.  In  fact,  the  data  set  could  have  been  reduced  by  one-half 
and  only  252  of  the  total  information  which  constrains  the  model  variableo 
would  have  been  lost. 

Transversely  polarized  shear  waves  (SH)  from  v;ell  located  earth¬ 
quakes  are  used  to  determine  an  upper  mantle  Ghcar  velocity  structure 
compatible  with  both  travel  tinea  and  waveforms  of  observed  seismograms. 
Vertically  polarized  shear  waves  (SV)  arc  shown  to  be  inadequate  in 
identifying  triplications  because  of  contamination  by  PL-couplcd  shear 
waves  and  other  P-SV  intcractiono.  The  observed  Si!  travel  times  arc 
considerably  slower  than  the  JB  values  out  to  30  degrees. 
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I  •  iunmary 


This  year's  research  In  Seismic  Phenomena  Connected  with  Earthquakes 
and  Explosions  has  concentrated  on  four  task  areas. 

1)  Evaluation  of  seismic  discriminants 

2)  Understanding  and  determination  of  seismic  source 

characteristics 

3)  Determination  of  the  earth's  seismic  properties  (i.e. 

velocity  and  density  structure  and  Q~^) 

4)  Inversion  of  observations  using  crust  and  upper  mantle 

properties  to  obtain  the  spatial  and  temporal  charac¬ 


teristics  of  earthquakes  and  explosions 
In  the  semiannual  technical  report,  we  reported  on  two  different 
approaches  for  investigating  the  characteristics  of  the  seismic  source,  in 
particular  the  determination  of  the  corner  frequency  for  P  and  S  waves.  The 
first  uses  theoretical  volume  source  models  of  tectonic  release  of  the 
Archambeau  type  and  obtains  numerically  their  P-  and  S-wave  spectra.  The 
spectra  for  these  particular  models  show  a  definite  difference  in  the  P-waves 
and  S-wave  corner  frequencies.  This  approach  yields  theoretical  Mg  vs 
curves  which  can  be  discussed  in  terms  of  scaling  by  the  various  source 
parameters.  The  second  approach  is  to  obtain  the  earthquake  P-wave  source 
history  from  comparison  of  theoretical  and  observed  seisn  .grams.  This  history 
is  then  used  as  the  S-wave  source  history  to  obtain  theoretical  seismograms 
which  are  compared  with  S-wave  observations  of  the  same  events. 

Two  possible  criticisms  of  the  first  approach  is  that  the  Archambeau  type  tec 
tonic  release  source  does  not  model  the  tectonic  region  realistically  and  that  the 
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presence  of  a  zero  frequency  minimum  on  the  far-field  spectrum  is  due  to  a 
questionable  termination  of  the  tectonic  release  region  at  a  finite  distance 
from  the  source.  Archambeau  claims  that  the  approximation  using  a  finite 
termination  is  realsistic  and  closer  to  actual  field  situations  than  using 
the  same  spatial  release  distribution  to  infinity.  We  are  continuing  research 
on  tectonic  release  sources  using  more  realistic  models  of  tectonic  stress 
distributions.  In  particular  we  are  trying  to  determine  for  what  stress 
a  sd  body  force  distributions  the  zero  frequency  minimum  is  a  valid  result. 

The  use  of  synthetic  body  and  surface  waves  is  becoming  an  extremely 
valuable  tool  in  not  only  determining  the  character  of  the  seismic  source 
but  also  in  defining  the  seismic  propagation  properties  of  the  earth.  In 
the  semiannual  report,  we  gave  an  example  of  a  method  using  amplitudes  and 
arrival  times  of  body  and  Rayleigh  waves  from  NTS  events  to  determine  an 
average  crustal  model  from  NTS  to  mid-Arizona.  In  this  report,  transversely 
polarized  shear  waves  (SH)  from  well  located  earthquakes  are  used  to  determine 
an  upper  mantle  shear  velocity  structure.  It  is  further  demonstrated  that 
observations  of  travel  times  for  vertically  polarized  shear  waves  (SV),  the 
most  common  technique,  can  lead  to  erroneous  shear  structure  models  dur  to 
continuation  of  P-SV  interactions. 

In  Section  III,  the  spatial  and  frequency  distribution  of  path-corrected 
Rayleigh  waves  from  the  San  Fernando  earthquake  are  systematically  related 
to  the  faulting  process  as  a  multiple  event.  The  surface  wave  source  is 
taken  to  be  a  depth-distributed  set  of  double  couples.  A  least-squares 
inversion  is  used  to  find  the  set  of  source  parameters  which  optimally  fit 
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the  variance-weighted  data.  In  this  study  the  strength,  the  slip  angle, 
and  the  rupture  delay  times  are  optimally  determined  for  the  spatially 
assumed  distribution  of  double  couples.  A  sophisticated  error  analysis  is 
performed  to  estimate  the  uncertainties  of  the  calculated  model  variables. 

One  of  the  interesting  side  results  from  this  inversion  technique  is 
the  a  priori  determination  of  what  azimuthal  distribution  of  observations 
is  most  beneficial  in  determining  the  source  parameters,  the  uncertainties 
in  the  calculated  model  parameters,  and  the  relative  coupling  between 
various  parameters.  Of  the  18  WWSSN  stations  used  for  this  event,  two 
stations  contain  a  total  of  302  of  the  total  information  in  the  data  set. 

In  fact,  the  data  set  could  have  been  reduced  by  one-half  and  only  252 
of  the  total  information  which-constrains  the  model  variables  would  have 
been  lost. 

One  could  visualize,  using  these  inversion  operators  to  design  multiple 
explosion  events  which  would  make  discrimination  criteria  unreliable  except 
in  azimuths  where  data  was  difficult  to  obtain.  Of  course,  the  inversion 
operators  used  here  were  for  Rayleigh  wave  data  only.  For  discrimination 
techniques  based  on  Love  and  Rayleigh  waves  possibly  combined  with  body  waves 
the  design  would  probably  be  extremely  complicated  and  economically  unfeasible. 

Of  equal  Importance,  this  section  describes  and  uses  techniques  for 
correcting  surface  wave  data  for  mixed  path  crustal  structure  between  source 
and  receiver,  and  for  multipath  propagation  interference.  The  error  bounds 
associated  with  these  corrections  are  estimated,  and  their  influence  on  the 
reliability  of  the  model  are  determined. 

Since  the  crust  and  upper-nantle  shear  structure  dominate  the  propa¬ 
gation  characteristics  of  Rayleigh  and  Love  surface  waves,  Section  IV  deals 
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with  an  attempt  at  constructing  a  preliminary  shear  velocity  model  compatible 
with  both  travel  times  and  waveforms  of  observed  seismograms.  Transversely 
polarized  (SH)  LRSM  and  WWSSN  observations  are  used  in  the  study.  Vertically 
polarized  shear  waves  (SV)  are  shown  to  be  inadequate  in  identifying  tripli¬ 
cations  because  of  contamination  by  PL-coupled  shear  waves  and  other  P-SV 
interactions.  Accurately  located  west  coast  earthquakes,  modeled  as  shear 
dislocations,  are  used  as  sources.  The  observed  SH  travel  times  are  con¬ 
siderably  slower  than  the  JB  values  out  to  30  degrees.  Select  profiles  of 
observed  waveforms  in  conjunction  with  the  travel-times  provide  the  data  for 
model  determinations.  This  is  accomplished  by  fitting  the  observed  waveforms 
with  synthetics. 

Research  is  continuing  on  the  relations  between  seismic  source  charac¬ 
teristics  and  magnitude.  Recently  Anderson  a.;J  Kanimori  have  made  a  study  on 
the  relationships  between  energy,  moment,  magnitude  and  source  dimensions. 
Using  both  deterministic  and  statistical  source  models,  they  have  shown  that 
the  empirical  Gutenberg-Richter  magnitude-energy  relationship  has  strong  theo¬ 
retical  support  and  implies  that  rupture  times  for  large  earthquakes  are  much 
longer  than  rise  times.  In  general,  earthquake  magnitudes  (M)  are  proportional 
to  L  ,  where  L  is  a  linear  source  dimension,  but  M  -  L3  for  small  events  and 
M  -  L  for  very  large  events.  Interplate  and  intraplate  earthquakes  form 
separate  populations  in  moment -magnitude,  or  moment-dimension  plots.  In  each 
population  the  inferred  stress  drop  and  effective  stress  is  roughly  constant. 
The  large  scatter  in  plots  involving  magnitude  is  considerably  reduced  when 
earthquakes  are  grouped  according  to  type.  The  scatter  virtually  disappears 

when  magnitude  is  removed  as  a  variable,  e.g.,  dimensions  versus  moment 
plots. 
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II.  Abstracts  of  Publications  and  Reports 
During  This  Contract  Period 
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Helmberger,  D.V.,  "Generalized  Ray  Theory  for  Shear  Dislocations,"  Bull. 
Seism.  Soc.  Am.,  in  press. 

Displacements  at  the  surface  of  a  layered  halfspace  produced  by  a  point 
source  dislocation  is  investigated.  Expressing  the  source  in  terms  of 
P,  SV  and  SH  displacement  potentials  allows  the  solution  to  be  expanded  in 
generalized  rays.  The  transient  response  for  each  ray  is  obtained  by  the 
Cagniard-deHoop  method.  First-motion  and  high-frequency  asymptotic 
approximations  of  the  exact  solutions  are  discussed. 

Numerical  experiments  on  the  far-field  response  produced  by  shallow 
dislocations  are  presented.  Synthetic  seismograms  computed  for  a  strike-slip 
dislocation  are  compared  with  some  observations  from  the  Borrego  Mountain 
earthquake.  Preliminary  results  indicate  an  equivalent  point  source  depth 
of  9  km  with  the  far-field  time  function  approximated  by  a  step  function 
with  an  exponential  decay.  This  time  function  fits  both  the  P  and  S  wave¬ 
forms.  The  apparent  shift  in  corner  frequency  between  the  P  and  S  waves 
for  shallow  events  as  reported  by  some  investigators  is  explained  by  surface 
reflections. 


Kurita,  T.,  "Upper  Mantle  Structure  in  the  Central  United  States  From  P 
and  S  Wave  Spectra1,"  Phys.  Earth  and  Planet.  Inter.,  in  press. 

In  order  to  study  regional  variations  in  the  upper  mantle  structure  in 
the  central  United  States,  t^e  second  step  of  a  two  step  procedure  formulated 
in  Kurita  (1972a)  has  been  taken  and  experimentally  shown  to  be  very  powerful 
method  for  elucidating  the  fine  configuration  of  the  low-velocity  zone. 

The  strong  advantage  of  this  method  as  compared  with  the  other  methods,  is 
a  mutual  independency  of  the  extent  of  the  velocity  decrease  and  the  depth 
to  the  bottom  of  the  low-velocity  zone.  The  upper  mantle  structure  down  to 
about  220  km  in  this  region  has  been  inferred  from  complex  transfer  ratios 
of  long  period  P  and  S  waves  from  deep-focus  earthquakes  recorded  at  FLO, 

OXF,  and  SHA  on  the  basis  of  the  crustal  models  obtained  in  Kurita  (1972b). 
From  the  Interior  Plain  to  the  Gulf  of  Mexico  the  low-velocity  zone  shifts 
to  a  shallower  dep'.n  while  increasing  its  thickness  and  decreasing  its 
velocity.  Thio  zone  is  about  50  km  thick  layer  ranging  in  depth  from  about 
150  to  200  km  under  the  Interior  Plain,  about  75  km  thick  layer  from  about 
120  to  195  km  under  the  Gulf  Coastal  Plain,  and  about  80  km  thick  layer 
from  about  95  to  175  km  under  the  continental  shelf  of  the  Gulf  of  Mexico, 
all  nearly  along  89°N  longitude.  The  decrease  in  S  wave  velocity  at  the  top 
of  this  zone  is  about  0.30,  0.45,  and  0.70  km/sec  under  each  of  the  areas 
above,  although  the  last  value  may  be  somewhat  an  overestimate.  Both 
boundaries  of  this  zone  are  sharp,  the  transition  occurring  over  at  most 
about  10  km.  In  this  region  the  existence  of  the  high-velocity  lid  zone 
is  possible. 
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Geller,  R.J.,  "Body  Force  Equivalents  for  Stress  Drop  Seismic  Sources," 
submitted  to  Geophys.  Journ. 

When  an  earthquake  point  source  is  represented  in  the  most  general 
way  (a  stress  drop  source)  six  independent  functions  of  time  (the  six 
stress  drop  elements)  are  required  to  describe  the  source.  The  common 
seismological  practice  of  using  a  single  time  function  for  a  point  source 
model  of  actual  earthquakes  is  not  valid  unless  all  six  time  functions 
are  independently  obtained  from  observed  data  and  found  to  be  the  same. 
The  difficulty  first  noted  by  Knopoff  and  Randall  (1970)  (that  one  double 
couple  lacked  sufficient  degrees  of  freedom  to  represent  an  arbitrary 
point  source  but  that  there  was  not  a  unique  decomposition  of  a  source 
into  a  combination  of  a  double  couple  and  a  compensated  linear  vector 
dipole  (CLVD))  is  shown  to  result  from  their  choice  of  non-orthogonal 
basis  elements.  The  unique  equivalent  force  representation  of  an  earth¬ 
quake  source  is  shown  to  be  three  perpendicular  force  doublets  parallel 
to  the  principal  axes  of  the  stress  drop  tensor.  The  success  of  the 
double  couple  source  model,  a  special  case  of  the  stress  drop  source,  in 
matching  most  observations  of  shallow  earthquakes  implies,  but  does  not 
require  that  most  shallow  earthquakes  result  from  shear  dislocations 
along  essentially  flat  surfaces. 


Geller,  R.J.,  "Representation  Theorems  for  an  Infinite  Shear  Fault," 
submitted  to  Geophys.  Journ. 

The  Green's  function  solutions  for  a  shear  dislocation  with  symmetry 
in  the  direction  of  infinite  length  have  been  integrated  analytically 
over  the  direction  of  infinite  length.  The  displacement  solution  is 
reduced  from  a  surface  integral  over  the  fault  plane  to  a  line  integral 
of  temporal  convolutions  over  the  width  of  the  fault.  If  the  time  history 
function  of  fault  displacement  is  any  piecewise  linear  function  the  con¬ 
volutions  can  be  integrated  analytically,  reducing  the  solutions  to  line 
integrals.  Numerical  results  and  plots  are  presented  for  a  simple  example. 
The  two-dimensional  solutions  give  exact  first  motions  at  points  "over" 
the  fault.  Also  the  solutions  can  be  used  to  test  two-dimensional  finite 
difference  and  finite  element  computer  programs. 


Helmberger,  D.V.  and  Engen,  G.R.,  "Upper-Mantle  Shear  Structure,"  submitted 
to  Journ.  of  Geophys.  Res. 

This  is  an  attempt  at  constructing  a  preliminary  shear  velocity 
model  compatible  with  both  travel  times  and  waveforms  of  observed  seismo¬ 
grams.  Transversely  polarized  LRSM  and  WWSS  observations  are  used  in  this 
study.  SV  signatures  are  demonstrated  to  be  inadequate  in  identifying 
triplications  because  of  contamination  by  PL-coupled  shear  waves  and  other 
P-SV  interactions.  Well  located  west  coast  earthquakes,  modeled  as  shear 
dislocations,  are  used  as  sources.  The  observed  SH  travel  times  are  con¬ 
siderably  slower  than  the  JB  values  out  to  30  degrees.  Select  profiles 
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of  observed  waveforms  in  conjunction  vi 
data  tor  model  determinations.  This  is 
observed  waveforms  with  synthetics. 


h  the  travel-times  provide  the 
accomplished  by  fitting  the 


«  J  “taTTSJ”  *  "T"  °f  models  as  preliminaty 

predict  strong  second  arrivals  beyond  30"  Modfllar8'’  ffrP  transltl°"s 
zones  between  400  and  700  km  nr^lL  w  M  d  ls  containing  low  velocity 
25".  None  o£  the  a£ove  features  C°”'plIcated  "aveshapes  between  20  and 

The  final  model  is  relatively  It  “PP31*"*  in  observed  SH  waveforms, 
as  the  HWB  P-model  except  that  the  ^  exhlbits  the  same  characteristics 
at  400  km  and  slightly  larger  at  SOn^6”™8**™10'1^  JUmp  ls  mUch  le8s 
similar  to  the  P-model.  m’  be  transition  is  subdued 


was'set'^ut  VrZ^TS  ^^Lr^Joenir^rL^na^0^;15^11  8el8B°8rflph  trailers 

ground  nuclear  test,  O.curo,  on  September  21  1079  r  Pre7an"ounced  under“ 

amplitude  information  were  used  tc^flna  .  1  Travcl  tl”e  and 

ing  synthetic  seismograms  using  the  CaenlardVHrn80  crus'al  l”odel  by  ralculat 
from  other  nuclear  events  at  NTS  ».  I  a  °P  rethod-  Sleigh  waves 

were  examined  as  a  con«ol  Le  mlnZ  r£  “  ^  T“CS°"  WSS  statid"' 

of  the  crust.  Croup  velocity  curves  we^8^  a  T'U"  °f  thC  top  halt 
waves  calculated  for  Tucson  a„S  hlngm'  tlism  ‘S.T  m'  Raylclgh 
teristics  of  Pn,  a  reflected  pn  -a  a  n  l  The  formations  and  charac- 

demonstrated  to'be  tolled of"^  r5l%  h"  :Xa"ined'  is 

and  contains  multiple  arrivals  of  p-cv  con  flectloa  from  the  mantle 

and  observed  seismograms  indicate  a  ctustal  £hicZss' Synthetlc 
Poisson  s  ratio  of  21  The  .  ctcness  or  JO  km  with  a 

jumping  from  6.1  to  7.9  km/sec.  Th^amoli appG3rs  to  be  sharP 
evidence  of  an;  lid  structure?  *  ^  3inplltude  behavior  of  Pn  shows  little 
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Hanks,  T.C.,  "Long-period  Strong  Ground  Motion  Following  the  San  Fernando 
Earthquake 

Two  hundred  thirty-four  components  of  ground  displacement  are  the  basis 
of  an  investigation  of  long-period  (approximately  2-10  seconds)  strong 
ground  motion  following  the  San  Fernando,  California,  earthquake  at  local 
(A  <  100  km)  distances.  Seismic  moment,  source  dimension,  radiation  pattern, 
source  propagation,  the  development  of  surface  waves  together  with  their 
subsequent  dispersion,  and  the  local  structure  of  the  Los  Angeles  Basin 
all  appear  to  have  first  order  significance  in  the  shaping  of  observed 
ground  displacements  and  their  variation  with  azimuth  and  epicentral  distance. 
Long-period  strong  ground  motion  occurred  with  considerably  greater  amplitude 
at  southern  azimuths  than  at  northern  azimuths;  at  comparable  epicentral 
distances  amplitude  differences  approached  a  factor  of  10.  Radiation  pattern 
alone  seems  insufficient  to  explain  the  observed  variations.  Especially 
strong  signals  at  azimuths  due  south  of  the  epicenter  and  the  development 
of  large  amplitude  (3-5  cm),  long  duration  surface  wave  trains  at  stations 
along  the  southeastern  margin  of  the  Los  Angeles  L^sin  suggest  the  signi¬ 
ficant  effect  of  this  geologic  structure  on  the  ooserved  ground  displacements. 
Data  redundancy  at  closely  spaced  stations  in  Los  Angeles  and  at  Lake 
Isabella  Dam  suggest  that  errors  in  ground  displacement  at  periods  near  8 
seconds  are  approximately  1  cm  and  at  periods  near  15  seconds  are  several 
cm,  for  a  record  sensitivity  of  7.6  cm/g. 


Newton,  C.A.,  "Temporal  Changes  in  P-Wave  Velocity  Near  Riverside,  California." 

P  travel  times  from  NiS  to  Goldstone  and  Riverside,  California  have  been 
n 

measured  for  the  last  eight  years.  During  this  period  there  have  been  four 
or  more  cycles  of  compressional  velocity  changes  in  the  crust  NNE  of  River¬ 
side,  but  not  necessarily  limited  to  that  area.  Periods  of  lower  velocity 
are  associated  with  decreased  seismicity  and  periods  of  higher  velocity  are 
associated  with  increased  seismicity.  These  data  suggest  that  this  area 
which  includes  portions  of  the  San  Jacinto  and  San  Andreas  fault  zones  have 
repeatedly  striined  to  a  dilatant  state  and  then  fractured.  The  largest 
earthquake  related  to  these  velocity  changes  was  a  magnitude  5. A  shock  which 
followed  a  low  velocity  period  of  approximately  450  days.  It  appears  that 
this  area  is  now  in  a  dilatant  state  which  began  on  August  1972. 

Whitcomb,  J.H.,  "P'P'  Times  and  Average  Mantle  Velocities  Under  Continents 
and  Oceans." 

P'P'  times  are  an  important  data  type  because  they  provide  one  of  the 
few  means  at  hand  for  comparison  of  mantle  structure  under  different  litho¬ 
spheric  regions  such  as  continents  and  oceans.  However,  scatter  of  the 
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data  for  this  phase  has  long  discouraged  investigators  from  using  the 
readings  for  anything  more  than  an  approximate  check  on  P'P'  times.  A 
partial  explanation  for  this  scatter  is  asymmetric  P'P'  which  arrives 
earlier  than  the  symmetrically  reflected  phat  s,  but  at  least  some  of 
the  difficulty  lies  in  assignment  of  the  incorrect  P'P'  branch  to  a 
reading  and  in  assumption  of  the  incorrect  velocity  structure  at  the 
reflection  point.  These  problems  are  avoided  and  the  analysis  suggests 
that  one-way  travel  of  P  waves  through  oceanic  mantle  is  delayed  by  0.65 
to  0.95  sec  relative  to  United  States  mid-continental  mantle. 


American  Geophysical  Union,  Fifth-Fifth  Annual  Meeting,  April  8-12,  1974, 
Washington,  D.C. 

Kanamori,  H.  and  Cipar,  J.J.,  "Focal  Process  of  the  1960  Chilean  Earthquake." 

Long-period  seismograms  recorded  at  Pasadena  are  used  to  determine  the 
focal  process  of  the  1960  Chilean  earthquake.  Synthetic  seismograms 
computed  for  various  fault  models  are  matched  with  the  observed  strain 
seismograms  to  determine  the  fault  parameters.  A  low-angle  thrust  model 
with  rupture  length  of  800  km  and  rupture  velocity  of  3.5  km/sec  is  consi¬ 
stent  with  the  observed  Rayleigh  to  Love  wave  ratio  and  the  radiation 
assymraetry.  A  seismic  moment  of  1  to  3  x  10  dyne  cm,  the  largest  ever 
reported,  is  obtained  over  a  period  range  from  200  to  SOO^sec.  This  value 
together  with  the  estimated  fault  area  of  about  2  x  10^km~  fits  into  the 
recently  established  moment  vs.  fault  area  systematics.  The  excitation 
efficiency  seems  to  increase  with  the  period.  The  strain  seismogram  clearly 
shown  unusually  long-period  (300  to  600  sec)  wave  arriving  at  the  P  time 
of  a  large  foreshock  whicr  occurred  15  minutes  before  the  main  shock, 
suggesting  a  large  slow  de  ormation  in  the  epicentral  area  prior  to  the 
major  failure;  the  entire  focal  process  may  be  envisaged  in  terms  of  a 
large-scale  deformation  which  started  rather  gradually  and  eventually 
triggered  he  foreshocks  and  the  "main"  shock.  This  mechanism  may  explain 
the  large  premonitory  deformations  decumented,  but  not  recorded  instru- 
mentally,  for  several  Japanese  earthquakes. 

Alewine,  R.W.  Ill,  "Application  of  a  Least-Squares  Inversion  Technique  to 
Surface  Wave  Spectra  to  Determine  Source  Parameters  of  a  Finite  Source." 

A  linear  inversion  scheme  is  presented  to  fit  in  a  least-squares  sense 
discretized  surface  wave  spectra.  With  this  technique,  an  estimate  of  the 
uncertainty  in  the  predicted  source  model  is  given  when  errors  in  the 
observed  data  are  considered.  An  earthquake  source  Is  represented  by  a 
series  of  double  couples  which  can  be  offset  both  spatially  and  temporally 
to  approximate  a  moving  rupture.  A  comparison  of  this  representation  of  a 
moving  source  and  that  classically  representing  a  finite  source  is  given. 

The  inversion  scheme  yields  the  source  parameters  (moment  and  slip  angle) 
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for  the  individual  double  couples  as  well  as  rupture  velocity.  Rayleigh 
wave  spectra  from  a  suite  of  stations  recording  the  1971  San  Fernando, 
California  earthquake  are  examined  and  the  resolvability  of  the  inverted 
source  parameters  are  discussed. 
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HI.  Least-Squares  Inversion  of  Surface  Above  Spectra  to  Obtain 
Source  Parameters:  The  1971  San  Fernando  Earthquake 
This  section  is  Chap*-  >r  5  and  Appendices  3  and  A  of  Ralph  Wilson 
Alewine's  thesis  (197A'l  titled  "Application  of  Linear  Inversion  Theory  Toward 
the  Estimation  of  Seismic  Source  Parameters".  References  and  details  cited 
in  this  section  can  be  found  in  the  thesis.  This  Thesis  has  been  submitted 
and  successfully  defended  as  partial  fulfillment  of  the  requirements  of  the 
Degree  of  Doctor  of  Philosophy  at  the  California  Institute  of  Technology, 
Pasadena,  California. 
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Chapter  5 

Least-Squares  Inversion  of  Surface  Wave  Spectra 
to  Obtain  Source  Parameters: 

The  1971  San  Fernando  Earthquake 

5.1  Introduction. 

In  this  chapter,  we  will  examine  the  spatial  and  fre¬ 
quency  distribution  of  Rayleigh  waves  from  the  I971  San 
Fernando,  California,  earthquake  and  how  this  distribution 
is  related  to  the  faulting  processes  of  thau  event.  The 
use  of  surface  waves  to  determine  gross  source  parameters 
has  been  a  common  tool  in  recent  seismological  investiga¬ 
tions,  e.g.  Wu  (1968),  Tsai  and  Aki  (1970,  1971),  Canitez 
and  Toksoz  (1971,  1972),  and  Mitchell  (1973).  Unfortunate¬ 
ly,  the  problem  of  systematically  finding  values  of  the 
source  variables  which  optimally  fit  the  observed  data  has 
not  yet  received  the  attention  due  the  problem.  With  the 
exception  of  a  Monte  Carlo  technique  by  Tsai  (1972)  and  an 
iterative  technique  by  Turnbull  et  al.  (1973),  the  "trial 
and  error"  method  was  usually  employed  in  finding  suitable 
variables  which  describe  the  data. 

We  hope  to  extend  the  usefulness  of  surface  observa¬ 
tions  by  showing  how  an  inversion  procedure  can  be  applied 
to  surface  wave  observations  in  order  to  see  what  informa¬ 
tion  can  and  cannot  be  gained  about  the  earthquake  source. 
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Al°ng  with  this  Inversion,  we  will  use  a  method  of  express 
ing  the  finiteness  of  the  fault  plane  by  numerically  sum¬ 
ming  point  source  solutions.  This  method  will  allow  for 
the  variation  of  parameters  over  the  fault  surface. 
Finally,  we  would  like  to  compare  the  solution  that  we  ob¬ 
tained  for  the  San  Fernando  source  mechanism  from  the  sta¬ 
tic  study  in  the  previous  chapter  to  that  obtained  by  the 
dynamical  surface  waves. 

5.2  Theoretical  Model. 

Spectra  from  a  Point  Source.  Theoretical  surface  wave 
spectra  will  be  computed  for  a  source  in  a  multilayered 
medium  which  approximates  the  San  Fernando  area.  The 
elastic  parameters  for  this  medium  are  given  in  Appendix  5. 
Far-field  amplitude  spectra  depend  on  the  surface  wave 
amplitude  response  of  the  layered  medium  and  also  the 
source  type,  depth  and  orientation.  Expressions  for  the 
frequency  dependent  media  response  factor,  which  is  depen¬ 
dent  of  the  source,  are  given  by  Harkrider  (1964).  Compu¬ 
tation  of  this  factor  involves  use  of  the  Thomson-Haskell 
matrices  for  the  multilayered  half-space  (Haskell,  1953). 
The  formulation  for  the  component  of  the  surface  wave 
spectra  which  is  due  to  different  source  types  and  the 
orientation  of  the  source  was  developed  by  Ben-Menahem  and 
Harkrider  (1964).  The  far-field  solution  for  Rayleigh  wave 
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spectra  for  a  point  double-couple  can  be  written  as 
S(w)  k  Ar  x(e,h)e"1(krT^T) 

A(w)  *  - j= -  ,  (5.1) 

Vkr 

where 

S(w)  *  transformed  spectral  source  function 
k  *  Rayleigh  wave  wave  number 
Ar  *  medium  response  factor 

r  =  radial  distance  to  receiver 

X(9,h)  *  complex  radiation  pattern  function,  where  e 

Is  the  azimuth  of  the  station  relative  to  the 
source  and  h  Is  the  source  depth. 

The  radiation  pattern  function  is  given  by  Ben-Menahem 
and  Harkrlder  (1964)  and  by  Harkrlder  (1970).  This  func¬ 
tion  Is  dependent  on  the  source  parameters  and  the  source 
depth.  At  a  given  depth,  this  azimuth  dependent  function 
Is  characterized  by  the  slip, or  rake,  angle  and  the  dip  of 
the  source-equivalent  fault.  The  numerical  techniques  used 
In  the  computer  programs  which  calculate  these  functions 
are  discussed  by  Harkrlder  (1970). 


The  source  function  is  assumed  to  be  a  step  function 
In  moment .  This  is  given  by 
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where  t  Is  the  origin  time  of  the  event. 

Source  Flnlteness  Effects.  In  the  above  discussion,  the 
source  was  represented  by  a  point  source.  However,  If  the 
dislocation  on  the  actual  fault  surface  initiates  at  one 
point  and  propagates  to  another  point,  the  fault  finiteness 
and  rupture  velocity  can  become  important.  Ben-Menahem 
(1961)  has  shown  how  the  dimensions  of  the  source  and  the 
speed  of  rupture  play  an  important  role  in  the  resulting 
spectral  radiation  patterns.  This  author  shows  that  spec¬ 
tra  from  a  source  which  radiates  energy  evenly  as  propaga¬ 
tion  occurs  was  modulated  by  the  shift  factor, 

2^  *1X  •  (5-3> 

where 


X  -  !cr  [§-  e°»«V3  . 

The  following  parameters  are  defined1 

w  *  angular  frequency 
C  «  Rayleigh  wave  phase  velocity 
b  =  horizontal  rupture  length 
v  ■  horizontal  rupture  velocity 

*  angle  measured  from  the  diiection  of  rupture. 
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However,  If  the  source  radiates  energy  unevenly  during 
p -opagation,  then  this  simple  modulation  no  longer  applies. 
One  such  method  of  uneven  radiation  can  occur  by  having  the 
dislocation  (i.e.,  moment)  vary  along  the  fault.  The  fact 
that  surface  -.rave  excitation  is  depth  dependent  means  that 
a  fault  which  has  a  vertical  component  of  propagation  can 
also  give  rise  to  uneven  radiation  during  propagation.  The 
analytic  modeling  of  such  rupture  would  be  very  difficult 
for  a  layered  medium,  however,  we  give  here  a  numerical 
approximation  to  this  phenomena. 

Numerical  Approximation  of  a  Propagating  Fault.  In  this 
section,  we  will  present  a  method  of  numerically  approxi¬ 
mating  the  surface  wave  spectra  due  to  a  propagating  rup¬ 
ture.  To  do  this,  we  will  consider  the  fault  to  be 
modeled  by  a  discrete  number  of  double  couple  sources 
arranged  spatially  to  reflect  the  length  of  a  fault.  Prop¬ 
agation  along  the  fault  occurs  by  "turning  on"  each  of  the 
sources  sequentially  at  a  time  given  by 


where  A  is  the  distance  to  the  ith  individual  source 
si 

from  some  origin  and  Vp  is  some  rupture,  or  propagation, 
velocity  that  we  wish  to  approximate.  The  spectra  from 
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each  of  the  sources  will  have  a  different  phase  at  some 
fixed  observation  point.  The  total  spectra  measured  at 
this  point  will  be  the  complex  sum  of  the  Individual  con¬ 
tributions.  The  phase  delay  of  the  Individual  propagation 
sources  relative  to  the  Initial  source  taken  to  occur  at 
the  origin  Is  given  by 

^1  **  R0 

T  =  — - -  +  T, 

P1  C(w)  dl 

where  RQ  Is  the  distance  from  the  Initial,  or  reference 

A 

source,  to  the  exterior  observation  point  and  R^  are  the 
distances  from  the  other  sources  to  the  observation  point. 
C (w)  Is  the  phase  velocity.  The  geometry  of  this  problem 
Is  illustrated  In  Figure  5.1.  From  this  figure,  we  can 
write 


"i  *  Cds"  +  Ro  -  2  aSiVos(*R-epn1/2, 

where  d  Is  the  horizontal  separation  of  the  Individual 
si 

sources  relative  to  the  initial  source  taken  In  the  direc¬ 
tion  of  rupture,  X^.  0p  is  the  "azimuth"  to  the  station 
relative  to  the  initial  source.  Appendix  4  goes  into  the 
details  of  how  the  spectra  are  summed  given  this  phase 
information.  Both  strike  slip  and  dip  slip  faulting  can 
be  approximated  with  this  geometry. 


Figure  5.1.  Geometry  of  the  assumed  multiple  double 
couple  source. 
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We  caution  here  that  the  angle  0p  as  defined  in 
Figure  5.1  Is  not  the  same  angle  as  the  angle  6  In  the 
complex  radiation  pattern  term  of  equation  (5.1).  A  pre¬ 
liminary  examination  of  Figure  1  In  the  Ben-Menahem  and 
Harkrider  (196*0  paper  would  indicate  that  these  two  angles 
are  equivalent,  out  their  definition  of  6  is  in  error.  The 
angle  0  should  be  measured  positively  clockwise  from  the 
positive  strike  direction  when  viewed  above  the  half-space. 

As  a  check  to  this  numerical  approximation,  we  will 
model  a  vertical  strike  slip  earthquake  and  check  the  re¬ 
sulting  spectra  with  that  calculated  by  the  method  of 
Ben-Menahem  (1961).  This  model  consisted  of  16  sources, 
each  at  a  depth  of  20  km,  with  dip  90°,  slip  angle  0°  (pure 
left  lateral  faulting),  and  identical  moments.  The  sources 
were  spaced  at  2  km  increments  along  the  strike  of  the 
origin  source  (Xs=l80°)  so  that  the  fault  length  and  prop¬ 
agation  length  were  30  km.  Spectra  were  computed  for  appa¬ 
rent  propagation  velocities  of  1.5,  3.0  and  *J.5  km/sec  at 
a  point  in  a  direction  *15°  to  the  direction  of  propagation 
in  order  to  avoid  a  node  in  the  radiation  pattern  of  strike 
slip  events.  These  spectra  were  then  compared  to  that 
calculated  for  a  single  point  source  with  the  shift  term 
for  this  fault  length  and  propagation  velocities.  The 
moment  of  the  single  source  was  taken  as  the  simple  sum  of 
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the  moments  In  the  discretized  model.  The  results  are 
shown  In  Figure  5.2. 

In  this  figure,  we  have  normalized  the  amplitudes  to 
the  values  calculated  for  a  non-propagating  point  source 
at  the  depth  of  the  source  at  the  point  of  Initial  rupture. 
By  utilizing  this  normalization,  we  cancel  the  media  source 
depth  effects  on  the  spectra.  The  ratio  shows  only  the 
effects  of  the  fault  finiteness.  In  Figure  5.2a,  we  see 
that  the  directivity  calculated  by  this  numerical  technique 
matches  very  closely  that  given  by  Ben-Menahem  (1961).  For 
propagation  velocities  3.0  km/sec  and  greater,  there  is 
less  than  1%  difference  in  the  two  calculations.  We  see 
that  the  directivity  is  a  very  pronounced  function  in  the 
period  range  of  interest  for  slow  propagation  velocities. 

Now  that  we  have  shown  that  the  directivity  function 
can  be  accurately  calculated  numerically,  we  now  examine 
the  effect  of  dip  slip  propagation  on  surface  wave  ampli¬ 
tude  spectra.  A  model  of  a  45°  pure  dip  slip  fault  was 
constructed.  The  initial  motion  occurs  at  a  depth  of  20 
km  and  propagates  at  a  constant  velocity  up  the  fault  to  a 
depth  of  2  km.  A  total  of  10  double  couple  sources  were 
spaced  at  a  2  km  depth  and  lateral  increments  to  approxi¬ 
mate  this  fault.  The  propagation  distance  along  the  fault 
is  28.3  km,  and  the  pure  horizontal  propagation  distance  is 


Period,  sec 

Figure  5.2.  Comparison  of  the  frequency-dependent 
directivity  function  given  by  the  Ben-Menahem  formu 
lation  and  the  numerical  approximation.  Upper — 
results  for  a  30  km  long  left  lateral  strike  slip 
fault.  Lower— results  for  a  28.3  km  wide  <*5°  dip 
slip  fault. 


Ill— 23 


-173- 

18  km.  The  spectra  were  calculated  at  an  exterior  point 
which  is  in  the  direction  of  propagation.  Spectra  were 
also  calculated  using  the  Ben-Menahem  formulism  for  the 
same  parameters.  All  the  amplitudes  were  normalized  to  a 
non-propagating  point  source  at  a  depth  of  20  km.  Figure 
5.2b  shows  a  comparison  of  the  resulting  directivities  for 
this  dip  slip  case. 

We  cee  in  this  figure  that  there  is  a  vast  difference 
in  the  spectra  computed  by  the  two  methods.  The  most  sig¬ 
nificant  difference  is  the  accentuation  of  the  amplitudes 
in  the  25-60  sec  period  range.  As  a  check  as  to  whether 
the  large  amplitudes  calculated  in  the  numerical  case  were 
due  to  a  Doppler  effect  or  due  to  the  non-uniform  depth 
excitation  along  the  propagation  length,  we  normalized 
these  amplitudes  to  the  spectrum  which  was  calculated  using 
the  numerical  approach  when  an  infinite  propagation  velo¬ 
city  was  used.  This  normalization  cancels  any  non-uniform 
depth  excitation  since  the  dtoth  distribution  was  identical 
in  both  cases.  The  resulting  directivity  was  almost  iden¬ 
tical  to  that  calculated  by  the  Ben-Menahem  method  and 
shown  in  the  lower  part  of  the  figure. 

We  thus  conclude  that  the  effect  of  a  source  radiating 
energy  as  it  propagates  from  one  depth  to  another  during 
faulting  can  be  much  larger  than  the  Doppler  effect  and 
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cannot  be  ignored  when  spectra  due  to  this  type  of  faulting 
are  being  modeled.  The  depth  effect  on  the  spectra  will 
vary  from  medium  to  medium  and  from  source  type  to  source 
type.  No  attempt  was  made  here  to  catalogue  this  effect 
for  various  media  and  propagation  parameters. 

Approximation  of  San  Fernando  Faulting.  In  the  previous 
chapter,  we  discussed  the  geophysical  observations  which 
constrained  the  fault  model  geometry  for  the  1971  San 
Fernando,  California,  earthquake.  In  this  chapter,  we 
would  like  to  adopt  an  approximation  to  this  same  geometry 
to  try  to  explain  the  radiation  of  surface  wave  energy  from 
this  fault.  This  approximation  is  accomplished  by  using 
four  double-couple  sources.  The  primary  constraint  for 
limiting  the  approximation  to  only  four  sources  was  econo¬ 
mic,  since  the  forward  problem  had  to  be  recalculated  after 
each  iteration  for  the  partial  derivative  matrix.  However, 
we  consider  that  the  approximation  chosen  is  adequate  in 
light  of  the  periods  of  the  surface  waves  used.  The  para¬ 
meters  which  describe  this  fault  geometry  and  remain  fixed 
in  the  inversion  are  listed  below. 


Source  No. 


1 

2 

3 

4 


depth(km) 

dip 

ds(km) 

14.5 

52° 

0.0 

10.5 

52° 

2.49 

6.25 

350 

7.51 

2.00 

350 

14.36 
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These  loci  of  the  four  point  sources  do  not  lie  on  a 
straight  depth-distance  line,  but  rather  define  the  center 
line  of  a  fault  plane  whose  dip  increases  with  depth.  In 
calculating  the  propagation  delay  time  for  each  of  the 
elements,  the  distance,  A^,  is  used.  A^  is  given  by 


(h 


-  h 


)2  +  d‘ 


These  distances  are  the  straight  line  distances  from  the 
initial  source  to  the  three  other  double-couples.  The 
delay  time  error  caused  by  using  this  distance  rather  than 
the  distance  measured  along  the  curved  "fault  plane"  is 
less  than  1  sec  in  the  extreme  case. 

5.3  Forward  Problem  Formulation. 

In  this  section,  we  wish  to  pose  a  relationship  be¬ 
tween  the  source  model  parameters  and  the  observed  spectr  i 
in  a  linear  form  as  in  equation  (2.3).  As  we  did  for  the 
static  problem,  we  will  propose  some  fault  geometry  .*nd 
use  an  inversion  scheme  to  solve  the  fault  model  parameters 
for  this  geometry.  In  order  to  describe  the  data,  we  will 
use  a  set  of  double-couple  sources  that  are  arranged  spa¬ 
tially  to  approximate  the  fault  model  system  which  was 
presented  in  the  previous  chapter.  Some  of  the  parameters 
of  the  individual  sources,  specifically  those  which 


III-26 


-176- 


describe  the  geometry  of  the  fault  system,  will  be  set 
initially  and  not  changed  in  the  inversion.  These  para¬ 
meters  are  the  depth  and  dip  associated  with  individual 
source  elements  and  the  spatial  relationship  between  ele¬ 
ments.  The  variable  parameters  for  each  element  are  taken 
to  be  the  moment  and  slip  angle.  In  addition,  we  will  in¬ 
clude  the  propagation  velocity  and  propagation  direction 
as  variables  for  the  entire  system. 

The  model  response  matrix,  A,  for  this  problem  is 
given  to  be  composed  of  partial  derivatives  of  the  form 


3  s, 

TX 


(5. M 


where  is  the  spectral  amplitude  value  at  some  frequency 
f^  (£*1,L)  measured  at  some  '’azimuth”  ep  (j*l,J)  with 

V 

respect  to  the  fault  system.  (k*l,K)  is  the  source 

parameter  which  is  allowed  to  be  varied.  The  matrix  A  has 
K  columns  and  L  x  J  rows,  where  K  is  the  number  of  model 
variables  (specified  below),  L  is  the  number  of  discrete 
spectral  points  at  a  single  station  and  J  is  the  number  of 
stations  used.  For  this  case,  L=10,  J*l8,  K=10. 

We  see  from  Appendix  4  that  the  addition  of  multiple 


sources  is  a  non-linear  operation.  This  problem  can  be 
overcome  by  mapping  only  small  Increments  of  the  mode]  Into 
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changes  In  the  data.  The  criterior  for  this  linearization 
Is  the  same  as  was  discussed  In  Chapter  2.  The  partial 
derivatives  which  compose  the  matrix  A  were  computed 
numerically  for  a  given  model.  This  was  done  by  the  simple 
differencing, 


m 


S^m+Am^)  -  S^m) 


(5-5) 


The  Incremental  values  of  the  source  parameters  used  were : 

?  6 

Mq  -  moment  «  0.1  x  10  dyne-cm 
AX  -  slip  angle  *  10° 

AXR  -  propagation  direction  *  10° 

AVr  -  propagation  velocity  ■  0.25  km/sec  . 


5 •  Inversion  Procedure. 

An  Inversion  procedure  will  be  adopted  to  find  the 
model  parameters  which  best  fit  the  observed  spectral  data. 
Since  the  number  of  data  points  far  exceed  the  number  of 
variables,  we  will  fit  the  data  In  a  least-squares  sense. 
Appendix  3  gives  the  derivation  of  the  Inversion  operators 
for  the  general  least-squares  case. 

Since  we  have  mixed  units  In  the  source  model  para¬ 
meters  and  thus  mixed  units  in  A  and  6m,  a  problem  of 
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scaling  arises.  One  way  of  correcting  for  the  dimension¬ 
ality  of  the  parameters  is  through  the  use  of  the  weight¬ 
ing,  or  model  correlation,  coefficients.  The  weights  of 
the  source  parameters  can  be  adjusted  so  that  the  weighted 
coefficients  of  the  rows  of  A  are  about  equal.  For  a  non¬ 
linear  problem,  this  scaling  would  be  model  dependent  and 
thus  would  probably  vary  in  the  inversion.  This  method, 
although  it  can  give  satisfactory  results  in  obtaining  a 
fit  to  the  data,  is  not  esthetically  pleasing.  A  satis¬ 
factory  method  was  found  to  solve  this  scaling  problem. 

This  was  accomplished  by  non-dimensionalizing  the  problem t 

A  6m  =  6d  , 

by  dividing  through  by  the  data  and  model  parameters .  The 
dimensional  matrices  and  vectors  are  replaced  by  their  non- 
dimensional  counterparts  in  the  inversion.  The  components 
of  the  quantities  are  defined  as: 


6d,  ■  —  6d, 

i  d  i 

°i 


6m,  ■*  —  6m,  . 

0  m ,  0 


We  note  that  this  scaling  is  also  non-linear. 
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Throughout  the  inversion  calculations  carried  out  in 
this  chapter,  we  will  limit  the  absolute  value  of  a  model 
perturbation  in  one  iteration  to  that  used  in  the  differ¬ 
encing  method  of  calculating  the  partial  derivatives  in  the 
forward  problem.  This  was  done  to  avoid  projecting  any 
model  perturbation  outside  the  range  of  assumed  linearity. 

As  a  test  to  the  inversion  scheme,  we  will  calculate 
the  theoretical  spectra  at  18  far-field  points  (correspond¬ 
ing  to  the  location  of  the  18  stations  which  are  described 
in  the  next  section)  for  a  given  set  of  source  parameters. 
The  non-variable  source  parameters  are  the  same  as  those 
used  to  describe  the  San  Fernando  faulting.  These  theo¬ 
retical  spectra  were  then  input  into  the  inversion  scheme 
as  data"  to  see  how  accurately  the  corresponding  variable 
source  parameters  could  be  recovered.  The  results  of  this 
test  art  shown  in  Figure  5.3.  The  initial  values  for  the 
source  parameters  are  that  shown  on  the  left  at  Iteration 
0.  The  source  values  used  in  calculating  the  "data"  are 
s^own  by  the  arrows  on  the  right  hand  side  of  the  figure. 

We  see  that  at  least  6  iterations  are  required  before  the 
solution  starts  to  converge  to  the  true  values,  and  after 
ll  iterations,  the  estimated  values  are  very  close  to  the 
true  values.  The  restriction  that  the  maximum  siie  of  the 
perturbations  can  be  no  larger  than  the  differencing 
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Figure  5.3.  Test  to  check  the  amplitude  spectra  in¬ 
version  technique.  Initial  values  of  the  fault  param 
eters  are  shown  on  the  left;  the  theoretical  values 
are  indicated  by  arrows  on  the  right. 
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values  slows  the  convergence  somewhat.  Not  shown  In  this 
figure  Is  the  estimate  of  the  propagation  velocity  and 
propagation  direction  for  each  iteration.  These  two  quan¬ 
tities  varied  more  smoothly  and  converged  more  quickly  to 
the  theoretical  values  than  did  the  moments  and  slip  angles 
The  propagation  velocity  and  propagation  direction  con¬ 
verged  to  *J.2  km/sec  and  120°  from  an  Initial  value  of 
3.2  km/sec  and  100°, respectively. 

In  this  test ,  the  variance  of  each  "datum"  was  taken 
to  be  5%  of  the  calculated  value.  The  test  was  repeated 
with  the  data  variances  set  at  25%  of  the  calculated 
values.  Although  the  convergence  to  the  true  model  values 
was  not  as  good  as  In  the  previous  example,  the  solution 
remained  stable  through  4'v,°  Iterations,  and  a  fair  estimate 
of  the  model  parameters  was  obtained.  The  results  from 
these  experiments  were  encouraging  enough  to  attempt  to 
apply  the  technique  to  actual  observed  surface  wave  spectra 
5 . 5  Data. 

The  source  of  the  surface  wave  data  used  In  this  study 
was  from  the  long  period  vertical  component  Instruments  of 
the  World-Wide  Standard  Seismograph  Network  (WWSSN).  Only 
the  vertical  Rayleigh  wave  component  at  each  station  was 
considered  In  this  study.  Data  from  the  Canadian  Standard 
Stations  and  the  WWSSN  stations  located  In  the  continental 
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Unit  ed  States  could  not  be  used  since  either  the  long 
period  vertical  component  was  Inoperative  or  off-scale 
during  the  passage  of  the  Rayleigh  wave  train.  The  original 
intent  in  collecting  this  data  was  to  limit  the  stations 
used  to  those  whose  epicentral  distances  were  less  than 
60°-70°.  However,  to  provide  a  more  uniform  and  extensive 
azimuthal  coverage  of  stations,  this  condition  was  relaxed 
somewhat.  The  closest  station  in  a  given  azimuthal  incre¬ 
ment  from  the  epicenter  which  gave  a  readable  record  was 
chosen.  A  total  of  18  records  was  chosen  to  be  examined, 
and  these  stations  provided  a  fairly  good  azimuthal  coverage 
for  this  event.  Table  5.1  gives  a  list  of  all  these  sta¬ 
tions  used,  the  epicentral  distances,  and  the  geodetic 
azimuth  of  the  station  with  respect  to  the  epicenter. 

The  Rayleigh  wave  signal  was  identified  on  each  of 
the  records  and  digitized  at  an  irregular  interval  taking 
adequate  care  to  sufficiently  define  the  tra^e.  The*-c  data 
were  then  linearly  interpolated  to  an  equal  increment  rate 
of  one  sample  every  two  seconds.  The  digitization  was 
restricted  to  a  velocity  window  of  between  about  4.3  km/sec 
to  about  1.5  km/sec.  There  was  a  slight  variance  to  this 
window  depending  on  the  distance  and  source-receiver  path. 
The  choice  of  the  velocity  window  was  based  on  a  visual 
Judgment  as  tc  when  the  signal  arrived  and  essentially 


TABLE 
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ended  at  a  particular  station  and  on  the  dispersion  curves 
expected  for  the  traveled  path.  Typical  of  the  dispersion 
curves  used  were  those  of  Brune  and  Dorman  (1963),  McEvilly 
(1964),  Santo  (1963),  Tryggvason  (1962),  and  Brune  (1969). 

The  Fourier  spectra  of  these  records  were  then  taken 
with  a  numerical  transform  routine.  The  spectra  were  cal¬ 
culated  for  the  frequency  range  from  60  sec  period  to  8 
sec  period,  with  512  discrete  frequency  points  In  this 
range.  The  spectra  were  smoothed  slightly  with  a  two-point 
Hanning  filter  such  that  for  the  1th  frequency  component  of 
a  particular  record 

fi(“>  -  |{jCf1.2(«)+f1+2(»)]+|[f1_1(u)+fl+i(B?]+fi(u)j  _ 

Instrument  Correction.  The  Instrument  frequency  and  ampli¬ 
tude  response  for  the  World-Wide  long  period  seismograph 
system  were  calculated  from  the  published  Instrument  and 
coupling  constants  for  that  Instrument.  These  calculations 
are  based  on  the  relations  given  by  Hagiwara  (1958)  and 
Mitchell  and  Landisman  (1969).  Both  an  amplitude  and  fre¬ 
quency  correction  were  supplied.  Since  some  of  the  sta¬ 
tions  that  were  used  in  this  study  operate  at  different 
instrument  gains  during  different  parts  of  the  year,  a 
check  was  made  on  each  seismogram  to  verify  the  listed 


Ill— 35 


-185- 


instrument  gain.  This  was  done  by  comparing  the  height  of 
the  calibration  pulse  to  the  applied  current  impulse.  The 
magnification  computed  in  this  fashion  was  within  about  lO# 
of  the  listed  magnification  except  for  station  SBA  which 
appeared  to  be  off  by  a  factor  of  2. 

Geometrical  Spreading  Correction.  All  of  the  amplitude 
spectra  A(w)  are  corrected  back  to  a  common  distance,  a 

0 

by  multiplying  by  the  well-known  relation  for  spreading  on 
the  surface  of  a  sphere 


'GS 


fi 


sin  A 


sin  A, 


(5.6) 


where  A  and  AQ  are  given  in  radians.  For  the  data  in  this 
case,  we  chose  our  reference  distance,  A  ,  to  be  2000  km. 
Attenuation  Correction.  Following  the  instrumental  and 
geometrical  spreading  corrections,  the  spectral  amplitudes 
were  corrected  for  attenuation.  This  correction  takes  the 
form 


Ac(u>)  =  A(w)  eY^  x 


(5.7) 


where  y(u>)  is  the  amplitude  attenuation  coefficient  appro¬ 
priate  for  the  particular  source  to  receiver  path,  and  x  is 
the  epicentral  distance.  Since  the  paths  from  the  San 
Fernando  epicenter  to  most  of  the  WWSSN  stations  used  in 


x  10“  km 


Figure  5.^.  Amplitude  attenuation  coefficients 
assumed  for  the  two  structures  used  in  this  study. 
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this  study  Involve  mixed  crustal  type  paths,  the  appropriate 
attenuation  coefficient  would  necessarily  be  a  path- 
averaged  value.  We  attempted  to  take  this  into  considera¬ 
tion  in  the  following  manner. 

We  assume  that  each  of  the  transmission  paths  can  be 
divided  into  t vo  types:  continental  and  oceanic.  Prom 
the  literature,  we  compiled  estimates  of  the  fundamental 
model  Rayleigh  wave  attenuation  coefficients  for  these  two 
media.  These  estimates  are  summarized  in  Figure  5.4.  The 
continental  attenuation  coefficients  are  reduced  from  data 
given  by  Gutenberg  (1945),  Gutenberg  and  Richter  (1936), 
Tryggvason  (1965),  Nuttli  (1973)  and  Mitchell  (1973) .  The 
oceanic  attenuation  coefficient  curve  was  reduced  from  data 
in  papers  by  Gutenberg  (1945),  Ben-Menahem  (1965),  and  Tsai 
and  Aki  (1969).  The  partitioning  of  continental  and 
oceanic  travel  paths  was  estimated  for  each  station,  and 
the  spectral  amplitudes  measured  at  that  station  were 
corrected  by  multiplying  by  the  following  relation, 

[xTyr (w)+xny  (w) ] 

C .  -  e  L  L  00  (5.8) 

A  y 

where 

YL(«)  =  the  continental,  or  land,  amplitude  dissipation 


coefficient 
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Yq(w)  *  the  oceanic  amplitude  dissipation  coefficient 
xT  ■  the  distance  traveled  over  a  continental  type 
path 

Xq  «  the  distance  traveled  over  an  oceanic  type  path. 

Table  5.1  lists  the  percentage  of  the  total  path  length 
that  is  spent  in  each  of  thg  assumed  attenuation  provinces. 
Also  listed  in  this  table  is  the  geometrical  spreading 
factor  given  by  the  equation  in  the  above  paragraph. 

Local  Crustal  Correction.  This  correction  attempts  to  com¬ 
pensate  for  the  effects  of  surface  waves  traveling  later¬ 
ally  from  one  structure  to  another.  In  the  case  under  con¬ 
sideration,  there  are  marked  variations  in  the  structure 
of  the  crust  and  upper  mantle  between  the  area  in  which  the 
surface  waves  were  generated  and  the  area  around  the  sta¬ 
tions  at  which  they  are  measured.  Fortunately,  the  correc¬ 
tions  applied  along  the  path  between  these  two  points 
cancel  out  so  that  only  a  local  crustal  correction  at  the 
receiver  needs  to  be  applied  to  correct  the  spectral  amp¬ 
litudes  back  to  the  medium  in  which  they  were  generated. 

An  approximation  which  has  given  satisfactory  results  in 
predicting  the  amplitude  changes  of  Rayleigh  waves  travel¬ 
ing  across  a  lateral  boundary  has  been  given  by  McGarr 
(1969)  and  McGarr  and  Alsop  (1967).  The  validity  of  this 
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approximation  has  been  verified  by  wave  propagation  studies 
in  laterally  heterogenous  media  using  the  finite-difference 
numerical  technique  (D.  Boore,  personal  communication, 
197*0. 

If  we  assume  that  the  total  energy  in  the  fundamental 
mode  Rayleigh  wave  remains  constant  in  the  process  of 
transmissions  across  the  boundary  margin,  then  we  can  write 


Je  .  Bsl"> 

As  V  - 

where 

Ag(w )  *  the  amplitude  in  the  source  medium 

Ar(u)  ■  the  amplitude  in  the  receiver  medium 

Ws(w)  =  the  normalized  energy  flux  in  the  source  medium 

WR(u)  ■  the  normalized  energy  flux  in  the  receiver 

medium. 

The  normalized  energy  flux  is  given  by 
W(w)  *  2  E(w)  U(w)  f 

where  E(w)  is  the  total  potential  or  kinetic  energy  at  a 
particular  frequency  normalized  to  the  surface  displace¬ 
ment.  The  calculation  of  this  energy  excited  in  a  given 
layered  earth  model  is  given  by  Harkrider  and  Anderson 


Energy  flux.  xIO  ergs/cm  -sec 
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5.5.  Energy  flux  as  a  function  of  frequency 

The'  Va^lous  crustal  and  upper-mantle  models. 
The  velocity-density  structures  for  the  various  models 
are  given  in  Appendix  5. 


* 


III-41 


-191- 

TABLE  5.2 


Crust  &  Upper  Mantle  Type 

Station 

San  Fernando  Crust  over 
Gutenberg  Continent 

- 

Source  Medium 

Greenlandic  Shield  with 
Ice  Cap 

- 

NOR,  KTG , 

SBA 

Iceland 

- 

AKU 

Irish  Continental  Shelf 

- 

VAL 

Iberian  Shield 

- 

MAL 

Mid-Atlantic  Ocean 

- 

BEC,  SJG 

Western  America  Tectonic 

- 

CAR,  BHP , 

COL 

Pacific  Ocean  East 

- 

GIE 

Pacific  Ocean  West 

- 

AFI ,  KIP, 

GUA 

Andean 

- 

PEL 

Central  Japan 

— 

MAT 

Table  5.2.  Crust  and  upper  mantle  structure  used 

for  the  local  crustal  structure  spectral 
amplitude  correction.  Velocity-density 
models  for  the  different  structures  are 
given  in  Appendix  5. 
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(1966,  equation  1).  U(<u)  is  the  group  velocity  for  the 

particular  medium. 

At  each  of  the  18  WWSSN  stations  used  in  this  study, 
the  crustal  structure  at  that  station  was  approximated  and 
the  velocity  dispersion  and  the  spectral  energy  density 
calculated.  The  crustal  models  used  in  this  calculation 
and  the  references  for  these  models  are  given  in  Appendix 
5.  The  energy  flux,  W(w),  given  by  these  calculations  is 
summarized  in  Figure  5*5*  The  stations  for  which  each 
crustal  model  correction  was  used  are  listed  in  Table  5.2. 

In  order  to  correct  the  spectral  amplitudes  back  to  a 
common  medium,  we  multiply  by  a  correction  which  normalizes 
the  amplitudes  to  the  source  medium.  This  correction  takes 
the  form 

(5.9) 

Multipath  Propagation  Removal.  The  complex  cepstrum  tech¬ 
nique  which  was  established  by  Schafer  (1969)  has  proven 
to  be  an  effective  means  of  removing  spectral  modulatioi s 
caused  by  the  interference  of  simultaneously  arriving 
signals.  Two  or  more  signals  can  arrive  simultaneously  at 
a  given  station  due  to  either  multipathing  of  a  single 
signal  or  can  result  from  the  addition  of  several  seismic 
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signals,  multisourcing,  or  some  combination  of  these  two 
effects.  In  this  study,  we  want  to  eliminate  as  best  we 
can  from  our  observed  spectra  that  modulation  which  is  due 
to  any  multipathing  effects  and  retain  any  modulations 
which  may  arise  from  the  addition  of  multiple  sources.  The 
spectra  which  have  been  de-modulated  from  multipathing 
effects  can  then  be  examined  by  this  method  to  try  to 
determine  the  existence  of  whether  the  source  is  a  single 
or  multiple  event. 

Cohen  (1970)  and  Flinn  et  al .  (1973)  have  used  the 
technique  to  detect  the  separation  of  the  body  phases  P  and 
pP  for  very  shallow  events,  while  Linville  (1971)  and  Tsai 
(1972)  have  applied  the  method  to  surface  waves.  In  the 
cepstrum  technique,  a  homomorphic  deconvolution  is  applied 
to  separate  the  components  of  a  convolution  of  a  seismic 
signal  and  a  multipath  operator.  Basically,  the  method 
acts  like  a  filter  applied  to  the  amplitude  spectra  where 
the  interferrence  effect  appears  as  a  scalloping  of  the 
spectra  over  a  wide  frequency  range.  When  a  long  pass 
filter  is  applied,  the  spectrum  modulating  effects  are  re¬ 
moved.  This  filtering  is  done  by  applying  the  appropriate 
bandpass  to  the  cepstrum.  Tsai  (1972)  and  Linville  (1971) 
give  a  summary  of  the  technique  involved  and  the  appro¬ 
priate  filter  bands  to  use  in  correcting  spectra  in  the 
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Figure  5.6.  Example  of  cepstral  filtering  at  two 

stations.  Original  spectra  are  shown  by  dashed  line* 

th?^ef5‘ned’  °r  flltered»  spectra  are  shown  by  the 
solid  line. 
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bandwidth  which  is  involved  in  this  study.  Such  a  long 
pass  filter  was  applied  to  the  observed  spectra.  This 
spectral  filtering  had  little  effect  except  to  smooth  the 
spectra.  A  typical  example  of  the  removal  of  the  spectral 
modulations  is  shown  in  Figure  5.6.  The  spectra  are  pre¬ 
sented  as  a  function  of  frequency  to  show  how  the  modula¬ 
tions  are  frequency  distributed.  In  this  figure,  the 
original  spectra  are  shown  by  the  dashed  line,  and  the 
spectra  after  the  application  of  the  cepstrum  technique 
are  shown  as  a  solid  line. 

We  saw  in  the  last  chapter  from  the  static  dislocation 
models  of  the  San  Fernando  earthquake  that  there  appeared 
to  be  two  areas  of  large  dislocation  —  one  in  the  hypo- 
central  region  and  one  near  the  surface.  If  the  disloca¬ 
tions  from  these  two  areas  acted  like  two  separate  sources, 
separated  spatially  and  temporally,  then  the  interx'erence 
of  these  two  sources  should  appear  in  the  complex  cepstrum. 
The  cepstra  were  calculated  at  several  of  the  stations  cho¬ 
sen  in  the  opposite  direction  of  the  apparent  rupture. 

This  was  done  to  maximize  the  effect  of  the  final  separate¬ 
ness  of  the  sources.  At  station  COL,  there  appeared  to  be 
some  signal  interference  at  a  delay  time  of  about  8  sec 
and  at  11.5  sec.  This  might  be  interpreted  as  the  inter¬ 
ference  of  the  signal  excited  at  the  hypocenter  by  another 
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signal  which  was  excited  when  the  fault  became  near  or 
broke  the  surface.  This  effect  was  not  as  apparent  at 
station  NOR,  so  we  feel  that  no  firm  conclusions  should  be 
drawn  from  this  particular  analysis. 

Corrected  Data.  All  of  the  above  corrections  were  applied 
to  the  measured  amplitude  spectra  at  the  WWSSN  stations. 
These  spectra  are  shown  in  Figure  5-7 •  The  upper  part  of 
this  figure  shows  the  location  and  azimuthal  distribution 
of  the  stations  used.  The  lower  part  of  the  figure  gives 
the  corrected  amplitude  spectra  at  each  of  the  stations. 

The  spectra  are  arranged  by  columns  with  geodetic  azimuth 
increasing  from  top  to  bottom.  The  spectral  values  indi¬ 
cated  by  a  dashed  line  in  the  figure  are  the  values  for 
which  the  total  correction  exceeded  a  factor  of  10.  These 
data  will  have  to  be  considered  less  reliable  than  the 
longer  period  data  because  of  the  uncertainties  in  the 
corrections.  Samples  of  the  spectral  values  were  chosen  at 
ten  discrete  periods:  60,  50,  40,  3^»  30,  26,  22,  20,  18 , 
and  16  sec.  It  was  felt  that  this  sampling  interval  was 
adequate  to  reflect  the  spectra  shape  at  the  individual 
stations.  Periods  shorter  than  16  sec  were  not  included 
into  the  inversion  data  set  because  of  the  larger  uncer¬ 
tainty  in  the  corrections  at  these  periods,  particularly 
in  the  attenuation  factor. 
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Estlmated  Errors  In  the  Dat-.a 


Errors  to  the  spectra  which 


we  have  measured  for  this  event  arise  from  two  sources: 
1)  errors  in  measuring  the  spectra  at  the  different  sta 
tlons  and  2)  errors  in  applying  the  corrections  to  the 


spectra. 

Errors  from  actually  measuring  the  spectra  could  arise 
rTOm  digitizing  errors  and  from  taking  too  short  a  signal 
to  get  all  the  spectral  information  at  certain  frequencies. 
Jne  station,  KIP,  was  redigitlzed  to  include  a  much  larger 
signal  and  the  spectral  amplitudes  recomputed.  The  recal¬ 
culated  amplitudes  were  found  to  be  at  most  3-H*  different 
on  the  average  than  that  originally  computed.  A  value  of 
3*  error  was  taken  as  the  estimate  for  this  error. 

By  far,  the  largest  error  involved  in  the  application 
of  the  spectral  corrections  arises  from  the  attenuation 


correction. 


If  we  assume  that  we  know  the  attenuation 


values  to  only  25*  of  the  values  given  in  the  figure,  then 
the  estimated  error  for  this  correction  would  be 


V“> 


We  next  can  assume  that  we  know  the  crustal  corrections  to 

only  25*.  This  turns  out  to  be  a  small  error,  except  at 

very  short  periods,  because  the  original  correction  is 
small. 
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Dlstrlbutlon  of  Data  Information.  In  Appendix  3,  we  have 
derived  an  operator,  J  ,  which  tells  how  the  information  in 
the  data  set  is  distributed.  The  use  of  this  operator  has 
been  demonstrated  by  Minster  et  al,.,  (1974) ,  and  these 
authors  give  a  lucid  description  of  its  properties.  The 
importance  of  a  particular  datum,  which  is  a  spectral  amp¬ 
litude  at  a  particular  frequency  measured  at  a  particular 
spatial  location,  depends  on  both  the  location  and  the 
accuracy  of  the  datum.  If  we  consider  each  WWSSN  station 
to  be  a  subset  of  the  entire  data  set,  then  the  relative 
importance  of  that  station  is  simply  the  sum  of  the  impor¬ 
tances  of  the  data  in  that  subset.  Figure  5-8  shows  how 
the  information  in  the  chosen  data  set  is  distributed.  In 
this  figure,  we  have  summed  the  relative  importances  of  the 
frequencies  at  the  individual  stations  and  plotted  this 
importance  as  a  function  of  the  geodetic  azimuth  of  those 
stations.  The  results  of  this  plot  are  somewhat  surpris¬ 
ing.  We  see  that  two  stations,  COL  and  KTG,  contain  a 
total  of  305?  of  the  total  information  in  the  data  set.  In 
fact,  the  data  set  could  be  reduced  by  half  and  only  25%  of 
the  total  information  which  constrains  the  model  variables 

would  be  lost . 

From  Appendix  3,  we  see  that  the  importance  of  a 
datum  does  not  depend  on  the  actual  value  of  the  datum 


Figure  5.8.  Distribution  of  data  information  with 
respect  to  summed  information  density  at  a  single 
station.  Information  density  is  normalized  such  that 
the  total  from  all  stations  is  100J. 
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lt^self  but  rather  the  estimated  error  In  the  datum  and  the 
model  response  operator,  A.  The  Importance  values  are 
dependent  on  the  source  model  only  In  the  sense  that  A  Is 
model  dependent.  Although  not  done  In  this  case,  this  type 
of  calculation  could  be  carried  out  before  the  data  set  Is 

gathered  to  determine  an  optimum  distribution  of  stations 
to  use. 

5.6  Inversion  Results. 

Fau!t  Model.  The  ten  model  parameters  which  describe  the 
source  system  were  found  by  the  least-squares  fit  to  the 
observed  spectral  data  measured  at  the  18  WWSSN  stations. 
The  best  fit  model  parameters  as  determined  from  the  in¬ 
version  procedure  are  as  follows: 


Source  H 

Depth 

(km) 

Moment 

(102^dyne-cm) 

Slip  Angle 
(^) 

1 

1^.5 

0.62 

66° 

2 

10.5 

0.20 

78° 

3 

6.25 

0.51 

82° 

2.0 

0.39 

70° 

propagation  direction  (\R)  =  990 
propagation  velocity  *  2.95  km/sec  . 


Figure  5.9  illustrates  the  degree  of  fit  to  the  data 
that  this  model  exhibits.  The  dashed  straight  line  in  the 
center  of  each  radiation  pattern  indicates  the  orientation 


[•] 


T=30  sec 


T=26  sec 


T=20  sec 


T= 16  sec 


F4cure  S  9.  Calculated  and  observed  radiation  patterns 
Orientation  of  the  fault  plane  (N70°W)  is  shown  by  the 
dashed  line  in  the  center  of  each  pattern.  xhe  amp.i. 
tudes  of  the  patterns  are  in  units  of  cn-sec. 
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of  the  strike  of  the  fault  system,  taken  to  be  in  the 
direction  of  N70°W.  Rupture  initiates  at  a  depth  of  14.5 
km  and  propagates  in  the  direction  of  S11°W.  The  fit  to 
the  long  period  data  (40-60  sec)  shows  that  the  radiation 
pattern  is  a  symmetric  dipole  pattern.  These  two  patterns 
show  very  little,  if  any,  signs  of  fault  finiteness.  The 
finiteness  of  the  fault  starts  to  become  apparent  in  the 
middle  two  figures.  For  these  periods  (26-30  sec),  the 
patterns  are  symmetric  in  a  direction  perpendicular  to  the 
direction  of  propagation  but  not  in  the  direction  of  prop¬ 
agation.  For  the  shorter  periods  -(16-20  sec),  the  effects 
of  the  rupture  length  and  non-uniform  excitation  become 
very  important  in  shaping  the  radiation  patterns. 

Uncertainty  in  Model.  We  would  now  like  to  estimate  the 
uncertainties  in  these  best  fit  model  parameters.  As  done 
in  Chapter  2  for  the  stochastic  inverse  case,  we  will  cal¬ 
culate  these  uncertainties  at  a  particular  confidence  level 
by  mapping  the  estimated  errors  in  the  data  into  errors  in 
the  model.  The  least-squares  solution  matrix  for  unweighted 
model  parameters  which  performs  this  mapping  is  the  matrix 
given  by  Mathews  and  Walker  (1964,  p.  366-367) 


Q  -  A*c  ■ _1A 

nn 


(5.10) 


We  note  here  that  this  solution  does  not  take  into  account 
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the  possible  weighting  of  the  different  model  parameters. 
This  is  no  problem  in  this  case,  or  in  any  other  case, 
since  the  problem  can  be  non-dimensionalized  as  we  have 
done  here.  The  non-dimensional  form  of  this  operator  has 
the  same  form  as  that  given  above. 

The  bilinear  product  of  the  matrix,  Q,  and  some  model 
perturbation,  q  ,  yields  a  confidence  region  of  resolvabil¬ 
ity.  That  is, 

q*Q  q  *  k2(c)  .  (5.11) 

This  equation  defines  the  hyperellipsoid  of  uncertainty  for 
a  particular  model  perturbation  at  a  confidence  interval, 
k(c) . 

By  fixing  k(c)  to  a  particular  value,  taken  to  be  1.96 
in  this  problem,  we  calculate  possible  model  perturbations 
for  which  the  equation  of  the  ellipsoid  is  solved.  Any 
model  perturbation  vector  which  lies  inside  this  ellipsoid 
is  unresolvable  at  the  confidence  interval,  k(c),  while  any 
model  perturbation  vector  which  protrudes  this  ellipsoid  is 
resolvable  at  that  confidence  level.  In  the  example  con¬ 
sidered  here,  we  will  limit  the  analysis  to  two  dimensions. 
With  this  limitation,  we  can  define  an  ellipse  of  uncer¬ 
tainty  for  combinations  of  model  parameter  pairs.  As  we 
will  see  in  the  examples  below,  this  analysis  also  allows 
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the  consideration  of  the  coupling  between  the  various 
model  parameters. 

The  ellipse  of  uncertainty  is  shown  for  various  fault 

model  parameter  pairs  in  Figure  5.10.  In  Figure  5.10,  the 

P  fi 

moments  are  given  in  units  of  10  dyne-cm.  and  the  angles 
are  expressed  in  degrees.  In  Figure  5.10a,  we  have  con¬ 
sidered  the  area  of  uncertainty  when  the  moments  of  the  two 

deepest  source  double  couples  are  allowed  to  vary.  We  see 

P 

that  M  has  an  uncertainty  of  only  .05  x  10  dyne-cm  when 
°1 

it  is  allowed  to  vary  alone.  Perturbations  larger  than 

this  value  would  be  resolvable  by  the  data  unless  the 

magnitude  of  the  moment  on  the  ad„ acent  double  couple 

(M  )  was  also  allowed  to  vary.  The  maximum  uncertainty 
°2 

that  Mq  can  take  in  this  coupled  system  is  shown  by  the 

projection  of  the  maximum  excursion  of  the  ellipse  onto  the 

M  axis.  The  maximum  uncertainty  of  M  is  similarly 
°1  °2 
calculated. 

The  ratio  of  the  uncertainty  measured  at  the  point 
where  the  ellipse  crosses  the  axis  of  a  model  variable  and 
the  maximum  uncertainty  of  that  variable  defines  a  coeffi¬ 
cient  which  reflects  the  coupling  of  the  two  variables. 

We  will  call  this  ratio  the  model  co-variance  coefficient 
(MCVC).  For  a  co-variance  coefficient  equal  to  1.0,  there 
is  no  coupling  between  parameters.  That  is,  the  maximum 


faSlt6mndp?  na”e)V  Uncertainty  ellipses  for  various 

batlon^whlch^al^lnslde'the  °f 
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figure  5*10  (f— i).  Uncertainty  ellipses  for  various 
fault  model  parameter  pairs.  Combinations  of  pertur¬ 
bations  which  fall  inside  the  ellipses  are  not  detect¬ 
able  by  the  data  at  the  95^  confidence  level.  Moments 
are  given  in  units  of  102o  dyne-cmj  angles  are 
expressed  in  degrees. 
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uncertainty  of  one  variable  is  not  dependent  on  the  uncer¬ 
tainty  of  the  other  variable.  For  coefficient  values  less 
than  1.0,  the  coupling  increases  with  decreasing  values. 
Figures  5.10a-c  show  how  the  coupling  between  the  indivi¬ 
dual  sources  decreases  with  Increasing  distance  between  the 
sources.  We  see  that  there  is  strong  coupling  between  M 

°1 

and  M  (MCVC*.30)  but  very  little  coupling  between  M 

2  0]_ 

and  M  (MCVO.85).  In  part  a  of  the  figure,  we  see  that 

the  coupling  is  in  a  positive  sense  between  M  and  M 

.  °1  °2 
l.e.,  an  increase  in  one  moment  is  traded-off  with  an  in¬ 
crease  in  the  other.  Between  M  and  M  ,  the  coupling  is 

°1  °3 

in  the  opposite  sense,  and  between  M  and  M  the  coupling 

°1  °4 

that  exists  is  in  a  positive  direction.  Figures  5.10d-e 
show  that  adjacent  slip  angles  are  slightly  coupled,  but 
much  less  so  than  their  respective  moment  values.  The  MCVC 
for  the  pair  is  .75.  As  the  distance  increases  be¬ 

tween  sources,  the  relative  coupling  between  the  slip 
angles  of  the  sources  decreases  rapidly.  Figures  5.10f-i 
show  how  the  moment  and  slip  angle  for  the  individual 
sources  are  related.  These  results  are  somewhat  surprising 
in  that  the  slip  angle  and  moment  for  individual  sources 
are  almost  completely  decoupled.  The  MCVCs  for  these  param¬ 
eter  pairs  range  from  .90  to  .95. 

In  these  figures,  we  observe  that  the  maximum  coupled 
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uncertainty  a  particular  variable  can  acquire  is  constant 
from  comparison  to  comparison.  This  uncertainty  is  the 
standard  error  for  that  variable.  Mathews  and  Walker 
(196^)  show  that  this  standard  error  is  given  by 


‘ii  ‘ 


Generalizing  this  expression,  we  can  say  that  the  values 
of  the  model  parameters  are  known  at  a  confidence  level  c 
when  the  errors  are  expressed  as 


k(c)  . 


(5.12) 


This  is  the  uncertainty  which  we  will  assume  for  the  model 
variables.  These  uncertainties  at  the  95?  confidence  level 
are  given  below: 


*  0.15  x  10  dyne-cm 

*  0.21  x  1026  dyne-cm 
=  0.13  x  102^  dyne-cm 
■  0.11  x  102^  dyne-cm 


18° 
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AX3  "  10° 

AX^  «  9° 

AXr  -  5° 

AVr  =  0.2  km/sec. 

Comparison  with  Static  Solution.  We  would  now  like  to  com¬ 
pare  the  solution  for  the  fault  model  of  the  San  Fernando 
earthquake  derived  with  surface  data  to  that  obtained  with 
the  static  displacements.  To  do  this,  we  will  sum  the 
moments  of  the  Individual  static  fault  elements  for  the 
elements  that  fall  within  a  depth  range  around  the  location 
of  the  double  couple  sources.  We  will  assume  that  the 
actual  static  moment  distribution  will  be  bracketed  by  the 
acceptable  dislocation  models  which  gave  the  maximum  and 
minimum  moments.  The  moments  fo"  the  two  models  are  shown 
In  Figure  5.11.  We  see  here  that  for  the  upper  10  km  of 
the  fault  system  there  Is  very  good  agreement  In  the  static 
and  dynamical  solutions.  The  only  significant  difference 
between  the  two  solutions  Is  the  moment  In  the  hypocentral 
region.  The  surface  wave  study  Indicates  a  moment  which  Is 
a  factor  of  2-3  times  greater  than  that  predicted  by  the 
static  data.  We  recall  from  the  previous  chapter  that  the 
uncertainties  in  this  area  were  larger  than  for  the  more 
shallow  areas.  In  the  surface  wave  study,  this  was  not  the 
case.  The  uncertainty  in  the  hypocentral  region  was  not 


Moment,  (xIO  dyne-cm) 
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Surface 

waves 


Figure  5.H.  Comparison  of  the  moment-depth  distri¬ 
bution  function  found  from  the  static  study  and  the 
surface  wave  study.  Static  moments  are  averaged  over 
the  indicated  depth  ranges.  The  static  moment  dis¬ 
tribution  from  the  dislocation  model  shown  in  Figure 
4.5  and  4. 1C  is  shown  as  a  solid  and  dashed  line 
respectively. 
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substantlally  different  from  the  other  areas.  If  we  assume 
that  the  fault  geometry  used  In  the  last  chapter  Is  correct 
then  the  surface  wave  study  Indicates  5.5  m  of  displacement 

p 

over  an  area  of  36  km  In  the  hypocentral  region.  This 
compares  favorably  with  that  derived  by  Hanks  (197M)  of 
7.5  m  over  about  the  same  area.  This  last  calculation  was 
based  on  still  another  data  set  of  this  event,  namely 
near-field  strong  motion  records. 

5.7  Conclusions. 

In  this  chapter,  we  have  seen  how  theoretical  spectral 
values  of  surface  waves  can  be  calculated  numerically  for 
an  approximation  to  a  finite  source.  We  have  shown  that 
when  non-uniform  excitation  of  energy  occurs  during  fault¬ 
ing  over  some  depth-distributed  fault  surface,  the  result¬ 
ing  spectral  amplitudes  can  be  far  different  from  that  of 
a  non-depth-distributcd  source  at  the  same  depth  of  ini¬ 
tiation. 

Rayleigh  wave  data  were  calculated  from  a  set  of  18 
WWSSN  stations  and  normalized  by  correcting  for  the  various 
source  to  receiver  paths.  An  inversion  scheme  was  used  to 
find  the  set  of  fault  model  parameters  which  fit  the  ob¬ 
served  spectral  data  in  a  least-squares  sense.  The  best 

p  c 

fit  solution  gave  a  distributed  moment  of  1.7  x  10  dyne- 
c;,i,  and  it  was  found  that  both  the  moment  and  slip  angle 
for  this  event  varied  along  the  down-dip  width  of  the  fault 
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The  slip  angle  varied  In  such  a  way  that  the  displacements 
became  more  predominantly  dip  slip  as  the  dislocation 
propagated  upward  from  the  point  of  Initial  rupture  at 
about  3.0  km/sec.  Also  Indicated  by  the  surface  wave  data 
Is  that  the  propagation  proceeded  In  a  direction  5°-10° 
away  from  the  normal  of  the  strike  of  the  fault.  The  fault 
model  solution  obtained  by  the  surface  wave  Inversion  com¬ 
pares  very  favorably  with  that  given  by  the  inversion  of 
the  static  displacement  data  with  the  exception  of  the 
initial  dislocation  in  the  hypocentral  region  where  the 
static  data  field  gives  poor  resolution. 

Wyss  (1971)  obtains  an  average  moment  from  the  long 

period  surface  waves  at  two  European  stations  (which  were 

not  included  in  this  data  set)  of  1.9  x  10  dyne-cm. 

Canitez  and  Toksoz  (1972)  fit  the  surface  wave  spectra  at 

6  stations  in  the  northeast  and  southeast  geodetic  azimuths 

with  a  single  point  source  at  a  depth  of  1^4  km.  These 

P  fii 

authors  obtain  a  moment  of  0.75  x  10  dyne-cm  and  a  rake 
angle  of  ^5° -  Obvious  reasons  for  obtaining  different 
values  of  moment  and  slip  angle  by  these  authors  and  the 
present  study  are  1)  the  non-depth-dlstributed  source  gives 
rise  to  a  different  radiation  pattern  than  a  depth- 
distributed  source,  and  2)  the  azimuthal  coverage  of  sta¬ 
tions  was  more  limited  in  the  Canitez-Toksoz  data  set  than 
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in  the  data  set  used  here.  Comparing  their  data  set  to 
that  used  here,  we  find  that,  If  the  spectra  at  the  sta¬ 
tions  used  in  those  authors'  study  were  discretized  in  the 
same  manner  as  in  this  chapter,  then  their  data  set  would 
have  only  approximately  40*  of  the  total  information  which 
constrains  the  model  variables  as  was  used  in  this  study. 

Finally,  we  feel  that  the  technique  of  the  hybridiza¬ 
tion  of  the  numerical  approximation  of  a  finite  fault  and 
the  least-squares  inversion  formalism  allows  the  maximum 
derivation  of  information  about  the  source  from  a  set  of 
corrected  surface  wave  observations.  Use  of  the  operators 
assoclaced  with  the  inversion  methods  provides  a  powerful 
tool  to  determine  a  priori  which  data  would  be  most  bene¬ 
ficial  in  determining  the  source  parameters  and  to  estimate 
the  uncertainties  in  the  calculated  model  parameters  and  to 
determine  the  relative  coupling  between  various  parameters. 
It  is  hoped  that  these  techniques  can  be  extended  to  the 
inversion  of  other  seismic  data.  The  simplest  extension 
would  be  the  inversion  of  Rayleigh/Love  spectral  ratios. 
This  type  of  data  would  yield  far  better  results  than  when 
each  is  analyzed  separately  because  of  the  cancellation  of 
error  inducing  path  effects. 
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Appendlx  3 

Derivation  of  Least-Squares  Inversion  Operators 

In  this  appendix,  we  will  derive  the  least-squares 
inverse  and  the  associated  operators  for  three  cases.  For 
convenience  with  respect  to  manuscript  preparation,  matrices 
will  be  given  here  in  regular  typed  form,  not  boldface  as 
in  the  rest  of  this  thesis.  This  practice  should  cause  no 
confusion  in  this  case. 

Case  I :  Least-Squares  Inverse  with  no  weighting  of  the 
model  components  and  no  inclusion  of  the  data 
variances . 

We  define  the  forward  problem  to  be 

A  6m  ■  fidQ  ,  (A3. 1) 

where  fid  is  the  observed  data  vector,  5m  is  the  true 
o 

model,  and  A  is  the  operator  which  maps  a  function  from 
the  model  space  to  the  data  space.  We  define  the  vector 
b^fim  to  be  the  best  estimate  of  the  model,  5m,  that  we  can 
obtain. 

In  this  case,  we  wish  to  minimize  the  fit  to  the  data 
such  that 


ej(b)  -  | | «d0  -  6d||2 


(A3. 2) 
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is  minimized.  Here  we  have  || Z | f 2  *  Z*Z,  and  6d  is  the 
r'a-L  ulated  fit  to  the  data  given  by 


Substituting  (A3. 3)  into  (A3. 2)  and  expanding,  we  get 

el(b)  *  «d0#<Sdo-6do#A  b-b«A»6do+b«A»A  b  .  (A3.H) 

Now  performing  a  first  order  perturbation  in  6b  of  (A3.4) 
and  setting  to  zero  implies 

b  "  ■  (''*A)+A*«do  ,  (A3. 5) 

which  is  the  result  we  are  seeking. 

Case_II:  Least -Squares  Inverse  with  no  weighting  of  model 
components  and  with  the  inclusion  of  the  data 
variance. 

Here  we  have  the  same  forward  problem  as  before,  however, 
each  datum,  dQ^ ,  has  associated  with  it  some  variance,  0\. 
If  we  assume  that  the  data  are  unbiased  and  are  statisti¬ 
cally  independent,  then  we  can  write  the  variances  in 
diagonal  form 
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We  want  to  minimize  the  fit  to  the  data  in  a  least  squares 
sense  such  that 

e?(b)  -  ||6d  -  6d| | |  <*3-6) 

1  u  nn 

is  minimized.  This  weighted  norm,  | | * | ! c  »  ls  deflned  fcy 

nn 

llzH§  'z'cn;lz- 

nn 

Substituting  (A3. 3)  into  (A3- 6)  and  expanding,  we  get 

cZ(b)-6d0*On;1SV6d0*Cn-1A  b-(Ab)«Cn;1{d0 
+  (Ab)*Cn'1(Ab) 

Performing  the  first  order  perturbation  of  6b  on  this  last 
equation  and  setting  the  result  to  zero,  we  obtain 

b  *  6m  =  ( AkCnn^A)^A  Cnnl6do  * 


(A3. 7) 
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Case  III;  Least-Squares  Inverse  with  weighting  of  the 

model  components  and  with  the  inclusion  of  the 
data  variance. 

For  some  problems,  it  might  be  useful  to  consider  a  set  of 
non-zero  positive  weights  for  the  model  components.  This 
weighting,  or  model  correlation  coefficients,  can  be 
defined  in  the  1  allowing  manner, 

\t1  .  .  .  .  0 


0  ’  ’  ’  •  WM 

The  error  induced  through  the  use  of  this  correlation 
coefficients  is  given  as  the  following, 

e2 (b )  *  b*Cs”1b  .  (A3. 8) 

The  error  of  the  fit  to  the  data  is  still  given  by 
equation  (A3. 6).  For  this  case,  we  want  to  minimize  some 
quadratic  sum  of  these  two  errors. 

e2(b)  *  e^(b)cos6  +  e2(b)sin6  . 
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2 

This  minimization  of  e  (b)  proceeds  as  before  giving 

b  =  6m  =  Q+(0)A#Cni"i16do  (A3. 9) 

where  Q+(0)  =  (A*C  ~1A  +  tan0C  “1)+ 

nn  s  s 

Substitution  of  (A3* 9)  into  (A3. 3)  enables  us  to  write 
6d  =  J  6dQ  . 

J  is  called  the  data  importance  matrix  and  is  given  by 
J  -  A  Qt(9)A,Cn'1  . 

The  data  importance  matrix  tells  how  the  information 
in  the  data  set  is  distributed  (Wiggins,  1972).  This 
operator  serves  the  same  role  in  the  least-squares  inverse 
as  does  the  model  response  operator  (equation  2.29)  in  the 
stochastic  inverse.  The  oata  importance  matrix  gives  us 
information  as  to  how  the  model  components  as  a  whole 
"see"  how  large  an  effect  individual  data  values  have  on 
the  estimation  of  a  model.  Minster  et  al .  (1974)  show 
that  the  trace  of  this  operator  is  Invariant  and  given  to 
be  the  number  of  independent  linear  combinations  of  the 
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data.  For  the  inversion  problem  to  be  determinant,  this 
sum  must  equal  the  number  of  independent  model  parameters. 

In  both  the  inversion  procedures  discussed  in  this 
thesis,  we  have  neglected  an  obvious  operator.  We  will 
call  this  operator  the  perturbation  relaxation  operator. 
These  operators  are  defined  below. 

6d  *  Pc  5d  (Stochastic  inverse) 

6m  «  PT  <,  6m  (Least-squares  inverse) 


where 


Pc  *  AWA* ( AWA*+C  tan9 ) + 
o  •  nn 

PL.S.  ’  <A‘CnnlA+Cs;ltane,+(ACnnlA) 

These  operators  snow  simply  that  in  the  stochastic  inverse 
case,  when  Cnn  has  large  values  for  certain  data,  these 
data  will  be  essentially  neglected  in  the  fit  for  the  sake 
of  fitting  other  better  determined  data.  In  the  least- 
squares  inverse  case,  when  C  has  a  very  small  value  for 
certain  model  components,  the  model  will  try  to  fit  the 
data  by  ignoring  the  contributions  from  these  particular 
model  components.  Also,  we  see  that  for  Cpn  and  Cgg  fixed. 
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the  point  along  the  trade-off  curve  influences  the  per¬ 
turbation  sizes  in  the  iterative  inversion  scheme. 
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Appendix  ^ 

Spectral  Summation  of  Two  Sources 
in  the  Frequency  Domain 

Consider  the  addition  of  two  seismic  sources  by 

spectral  summation.  The  first  source  is  represented  by  an 

i<J>l  (w) 

amplitude  and  a  phase  spectrum  of  the  form  A^(w)  e 

The  phase  spectra  4>^(w)  contains  information  about  both 

the  source  and  the  propagation  path.  To  this  spectra  we 

wish  to  add  an  additional  source  represented  by 
i[<t>P(u)  +  A<t>2(w)] 

Ap(w)  e  ,  where  A2(w)  is  the  amplitude 

spectra  of  the  secondary  source,  2 ( uj )  is  the  phase  spectra 
of  the  secondary  source,  and  A<j>2(w)  is  the  additional  phase 
of  the  secondary  source  caused  by  a  spatial  or  temporal 
offset  of  the  second  source  from  the  first  source.  Propa¬ 
gation  effects  for  traveling  a  slightly  different  path 
from  origin  to  receiver  are  taken  into  account  with  this 
term.  If  we  let  z  be  the  total  of  the  two  sources,  then 
we  can  write, 

i*U)  i<Mu>) 

z  =  7(w)  e  =  A.  (cj)  e  + 

1  (A4.1) 

i[<J>2(u)  +  A4>2 ( uj )  ] 


A2(w)  e 
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Thus  , 

1/2 

Z(w)  =  {aJ(u)  +  A2(u)  +  2A1(w)A2(u)cos[^1(u)-(>2(w)-A^2(w)]} 

( AH .2 ) 

and 

A1(h))sln[»1  (u»)]-*-A?(m)sln[(»2(u))^A»2(M)  ] 

4»(“)  *  tan  1  A1(oj)cos[4>1(aj)]+A2(u))cos[*2(u))+A<))2(u)] 

(AH.3) 

The  phase  delay  of  the  single  original  source  is 
given  by 

T  .  h(ll  (AH.lt) 

Pi  10 

and  the  new  phase  delay  of  the  summed  signal  is 

T  =  tial  .  (AH. 5) 

Pz  u 

The  group  delay  of  the  summed  signal  is  given  by 

T  =  iMifjj  _  (AH. 6) 

3  to 

In  practice  in  order  to  avoid  as  much  as  possible  the 
difficulties  encountered  by  the  multivaluedness  of  the 
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trigonometric  functions,  we  wish  to  compute  only  the  addi' 
tional  group  delay  the  additional  source  has  over  the 
group  delay  of  the  single  original  source.  If  we  define 
this  additional  group  delay  by  6t  ,  then 

e 


6t 


g 


Tg„  "  T 


(All. 7) 


where  t  is  the  group  delay  of  the  original  signal.  By 
1  3<f»2(u))  3A<p0(u>) 

deflnlng  -17“  *  TE2  and  — •  *g  .  -e  have  by 


equation  (All.  6) 


V- 


Al(w)xg  +A2(w)  2~  +  r+A1(“)A2(u))cos(a) 

T_  =  - - i— - 


8, 


where 


A^(u))+A2(aj)+2A^(a))A2(w)cos(a) 


(Ail. 8) 


a  -  <P2  +  A<f.2  -  <t>1 

r+  =  T„  +  T  +  T 

+  El  S' 


g2s 


r  =  Tcr  +  T 

Ep  g 


-  T. 


2s  gl 


Substituting  this  expression  into  equation  (A4.7)  we  get, 
2 

A^Cw)  r_  +  A^(w) Agtw)  cos(a)  r 

6t  =  — - - - - - - 

A^(a))  +  A2(w)  +  2A^(uj)A2^*<)cos(a) 


(All. 9) 
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The  additional  phase  delay  caused  by  the  spatial  and 
temporal  offset  of  the  secondary  source  is  given  by 


T 


P 


2s 


t  T 

<*>  cWT  d2  • 


(A4.10) 


Here  R  is  the  distance  from  the  hypocenter  of  the  first 

signal  to  the  observer  and  R  is  the  distance  from  the  hypo- 

center  of  the  second  signal  to  the  observer.  C(w)  is  the 

phase  velocity  of  the  medium.  T,  is  the  delay  in  time 

a2 

after  the  initiation  of  the  first  source  before  the  origin 

of  the  second  source.  This  time  delay  can  be  expressed  as 

R 

an  apparent  rupture  velocity  given  by  T  ,  =  ~  f  where  R 

d2  VR  s 

is  the  separation  distance  for  the  two  events. 

Differentiation  of  (x  •  w)  with  respect  to  w  yields 

p2s 

the  secondary  group  delay,  x  given  by 

'  S2s 

A 

xg2s  ‘  cfSy  +  \  •  (ah. id 


where  U(w)  is  the  group  velocity  of  the  medium. 

It  is  seen  that  calculation  of  the  group  delay  does 
not  ii  /olve  any  Inverse  trigonometric  functions  so  that  the 
multivaluedness  of  the  functions  is  not  important.  However, 
this  is  not  the  case  for  the  calculation  of  i|/(w).  Jn 
practice,  the  values  of  ip(u)  are  numerically  differenced 


to  obtain  an  approximate  expression  for  the  group  delay. 

If  this  approximate  group  delay  is  more  than  one  period 

different  from  the  exact  calculated  group  delay,  then  a 

period  is  added  to  the  phase  delay,  x  ,  and  the  approxi- 

pz 

mate  group  delay  is  recalculated.  By  starting  with  the 
phase  and  group  delays  of  the  long  period  end  of  the  spec¬ 
trum,  the  true  phase  can  be  unwound  from  that  calculated 
numerically . 

For  the  addition  of  many  different  sources,  each  off¬ 
set  in  time  and  space  with  respect  to  one  another,  the 
above  calculations  are  repeated  for  each  source,  with  the 
original  amplitude  and  phase  spectrum  replaced  by  the 
partial  sum  spectra  after  each  individual  summation. 
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Appendix  5 

Velocity-Density  Models  for  Local  Crustal  Corrections 

This  appendix  lists  the  multi-layered  velocity-density 
models  of  the  crust  and  upper  mantle  used  in  the  computa¬ 
tion  of  the  local  crust  il  corrections  for  the  WWSSN  sta¬ 
tions  used  in  the  surface  wave  study  in  Chapter  5. 

The  following  symbols  are  employed  in  the  tables: 

D  =  layer  thickness  in  km 

Alpha  =  compressional  velocity  in  km/sec 

Beta  =  shear  velocity  in  km/sec 

Rho  =  density  in  gram/cm3. 

References  for  the  origin  of  the  data  used  in  the 
models  are  given  below  each  table. 
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San  Fernando  Crust  over  Gutenberg  Continent 


D 

Alpha 

Beta 

Rho 

0.50 

2.50 

1.20 

2.10 

3.00 

3.80 

2.50 

2.50 

2.00 

5.50 

3.20 

2 . 60 

14.50 

6.06 

3.40 

3.00 

15.00 

6.70 

3.80 

3.00 

13.00 

7.96 

4.60 

3.37 

25.00 

7.85 

4.50 

3.39 

50.00 

7.37 

4.21 

3.40 

75.00 

S.00 

4.41 

3.45 

50.00 

8.20 

4.50 

3.47 

100.00 

8.40 

4.60 

3.50 

100.00 

9.00 

4.95 

3.63 

100.00 

9.63 

5.31 

3.89 

Source:  Jungels  (1973) >  Gutenberg  (1944),  Wood  and 
Richter  (1933),  Harkrider  et  al .  (1963). 
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Greenlandic  Shield  with  Ice  Cap 


D 

Alpha 

Beta 

Rho 

2.50 

3.93 

1.9*4 

0.91 

16.50 

6.25 

3.7*4 

2.80 

23.70 

6.60 

3.85 

2.85 

37.30 

8.05 

*4.67 

3.30 

25.00 

8.10 

*4.72 

3.30 

100.00 

8.20 

*4.5*4 

3.  *4*4 

100.00 

8.30 

*4.51 

3.53 

80.00 

8.70 

*4.76 

3.60 

100.00 

9.30 

5.12 

3.76 

Gregersen 

(1970), 

Brune  and  Dorman 

(1963) 

Iceland 


D 

Alpha 

Beta 

Rho 

4.50 

4.70 

2.70 

2.60 

3.50 

•6.30 

3.60 

2.80 

10.00 

6.71 

3.90 

3.08 

10.00 

7.38 

4.30 

3.15 

50.00 

8.00 

*4.68 

3.30 

160.00 

7.90 

14.140 

3.35 

100.00 

8.20 

4.58 

3.1*0 

200.00 

8.20 

14.68 

3. HO 

Source:  Tryggvasor.  (1962),  B§th  (i960),  Tyggvason  and 
B&th  (1961),  Tryggvason  (1973). 
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Trish  Continental  Shelf 


D 

Alpha 

Beta 

f  .50 

1.52 

0.00 

1.50 

2.20 

1.57 

2.00 

5.40 

3.18 

20.00 

6.10 

3.55 

6.00 

7.30 

4.21 

50.00 

8.10 

4.61 

100.00 

8.10 

4 . 40 

100.00 

8.10 

4 . 60 

200.00 

8.10 

4.89 

:  Blundell  and  Parks  (1969),  Bamford 


Rho 

1.03 

2.55 
2.70 
2.82 
3.10 
3.30 
3.40 

3.55 

3.60 

(1971). 
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D 

Iberian 

Alpha 

Shield 

Beta 

Rho 

2.00 

3. 40 

2.00 

2.30 

18.00 

5.90 

3.50 

2.80 

10.00 

6 . 60 

3.70 

2.90 

20.00 

7.60 

4.50 

3.30 

30.00 

8.10 

4.70 

3.35 

100.00 

8.15 

4.20 

3.40 

100.00 

8.49 

4.77 

3.53 

100.00 

8.81 

4.89 

3.60 

200.00 

8.81 

4.89 

3.60 

Source:  Payo  '1970),  Payo  (1964). 
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Mid-Atlantic  Ocean 


D 

Alpha 

Beta 

Rho 

4 . 50 

1.51 

0.00 

1.03 

0.05 

1.52 

0.15 

1.65 

0.10 

1.60 

0.19 

1.70 

0.10 

1.71 

0.37 

1.79 

0.10 

1.80 

0.53 

1.36 

2.00 

5.00 

2.85 

2.66 

4.50 

6.69 

3.90 

3.06 

50.00 

3.00 

4.68 

3.30 

160.00 

7.90 

4.60 

3.35 

200.00 

8.20 

4.68 

3.40 

Source:  Ewing  (1969),  Talwani  et  al.  (1965),  Ewing  et  al. 
(1966),  Katz  and  Ewing  (1956),  Officer  et  al. 
(1952). 
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Western  America  Tectonic 


D 

Alpha 

Beta 

2.50 

3.00 

1.73 

24.50 

6.20 

3.5*1 

13.00 

6. 80 

3.87 

45.00 

7.80 

4.25 

50.00 

8.20 

4.38 

100.00 

8.20 

4.38 

100.00 

8.20 

4.50 

100  00 

8.20 

4.70 

:  Kanamori 

(1970), 

Bucher  and  Smith 

Alexander 

(1963), 

Ewing  and  Press 

Heezen  (1972) . 


Rho 

2.40 

2.83 

2.99 

3.30 

3.43 

3.52 

3.57 

3.62 

(1971), 

(1959), 
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Pacific  Ocean  East 


D 

Alpha 

Beta 

Rho 

3.80 

1.5i 

0.00 

1.03 

0.01 

1.52 

0.15 

1.65 

0.10 

1.60 

0.19 

1.70 

0.10 

1.71 

0.'? 

1.79 

0.10 

I.80 

0.53 

1.36 

1.30 

5.07 

2.8G 

2.67 

5.00 

7.20 

3.98 

2.90 

50.00 

8.00 

4.68 

3.30 

60.00 

7.90 

4.60 

3.35 

100.00 

7.90 

4.60 

3.35 

200.00 

8.20 

4.68 

3.40 

Source:  Piernattei  and  Nowroozi  (1969),  Santo  (1963), 

Saito  and  Takeuchi  (1966),  Dorman  et  al .  (I960), 
Raitt  (1956),  Ewing  et  al ■  (1969),  Latham  and 


Nowroozi  (1968),  Kovach  and  Press  (196I). 
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Pacific  Ocean  West 


D 

Alpha 

Beta 

Rho 

5.50 

1.51 

0.00 

1.03 

0.10 

1.60 

0.30 

1.65 

0.70 

5.07 

2.88 

2.67 

4.50 

7.20 

3.98 

2.90 

50.00 

8.20 

4.72 

3.30 

60.0C 

7.90 

4.40 

3.35 

100.00 

7.90 

4.40 

3.35 

o 

o 

o 

c 

C  J 

8.50 

4.74 

3.40 

:  Piermattei  and  Nowroozi  (1969),  Santo  (1963), 
Seito  and  Takeuchi  (1966),  Dorman  et  al .  (I960), 
Raitt  (1956),  Ewinf  et  al.  (1969),  Latham  and 
Nowroozi  (1968). 


Source 
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Andean 


D 

Alpha 

Beta 

Rho 

5.00 

5.00 

2.88 

2.67 

15.00 

6.00 

3.46 

2.87 

25.00 

6.6  0 

3.58 

2.99 

100.00 

7.90 

4.45 

3.30 

100.00 

8. CO 

4.40 

3.32 

100.00 

8.00 

4.4o 

3.32 

100.00 

8.20 

4.67 

3.38 

100.00 

8.20 

4.68 

3.38 

:  James  (1971),  Fisher  and  Raltt  (1962),  Clsternas 

(1961). 
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Central  Japan 


D 

Alpna 

Beta 

Rho 

1.00 

2.50 

1.47 

2.50 

5.00 

5.50 

3.18 

2.70 

16.00 

6.00 

3.40 

2 . 80 

10.00 

6.50 

3.70 

3.00 

30.00 

7.70 

4.37 

3.20 

100.00 

8.00 

4.50 

3.30 

200.00 

8.00 

4.50 

3.30 

Source:  Kaminuir.a  (1966),  Mikumo  (1966),  Asada  and  Asano 
(1972),  Kurlta  (1971). 
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IV.  UPPER-MANTLE  SHEAR  STRUCTURE 
Donald  V.  Helmberger  and  Gladys  R.  Engen 
Set  sinolog  leal  Laboratory 
California  Institute  of  Technology 
Pasadena,  California  91109 

Abstract 

This  is  an  attempt  at  constructing  a  preliminary  shear  velocity 
model  compatible  with  both  travel  times  and  waveforms  of  observed  seismograms. 
Transversely  polarized  LRSM  and  WWSS  observations  are  used  in  this  study. 

SV  signatures  are  demonstrated  to  be  inadequate  in  identifying  triplications 
because  of  contamination  by  PL-coupled  shear  waves  and  other  P-SV  interactions. 
Well  located  west  coast  earthquakes,  modeled  as  shear  dislocations,  are  used 
as  sources.  The  observed  SH  travel  times  are  considerably  slower  than  the 
B  values  out  to  30  degrees.  Select  profiles  of  observed  waveforms  in  conjunction 
with  the  travel-times  provide  the  data  for  model  determinations.  This  is 
accomplished  by  fitting  the  observed  waveforms  with  synthetics. 

Synthetics  are  computed  for  a  number  of  current  models  as 
preliminary  attempts  at  fitting  the  data.  Models  containing  large,  sharp 
transitions  predict  strong  second  arrivals  beyond  30°.  Models  containing 
low  velocity  zones  between  400  and  700  km  predict  complicated  waveshapes 
between  20  and  25*.  None  of  the  above  features  are  apparent  in  observed 
SH  waveforms.  The  final  model  la  relatively  smooth  and  exhibits  the  same 
characteristics  as  the  HWB  P-model,  except  that  the  percentage  velocity 
lump  is  much  leas  at  400  km  and  slightly  larger  at  500  km.  The  600  km 
transition  is  subdued  similar  to  the  P-model. 

^Contribution  No.  2460  ,  Division  of  Geological  and  Planetary  Sciences, 

California  Institute  of  Technology,  Pasadena,  California  91109. 


INTRODUCTION 


Considerable  effort  In  recent  years  has  gone  into  Interpreting 
upper  mantle  structure  in  terns  of  mlneraloglc  and  chemical  conposltions. 

A  fundamental  contribution  to  this  subject  was  given  by  Anderson  [1967) 
who  suggested  that  the  smooth  transition  over  the  413  to  1000  km  depth 
interval  described  by  the  Jeffreys  model  was  incorrect  and,  that  Instead, 
t:‘e  velocity  increased  rapidly  over  two  regions  near  depths  of  400  and  600  km. 
He  interpreted  these  transitions  in  terms  of  an  ollvlne-to-spinel  phase 
change  at  400  km  and  a  splnel-to-post  spinel  near  600  km. 

More  recently,  McGetchln  and  Silver  (1970]  and  others  demonstrated 
the  likelihood  of  other  coexisting  minerals  such  as  garnet  and  pryoxene. 

The  rather  complicated  problem  of  predicting  the  changes  of  seismic  properties 
across  phase  boundaries  involving  a  multi-mineral  assemblage  has  been  attempted 
by  Ahrens  [1973],  By  varying  the  temperature  and  the  iron  content,  assuming 
•i  uniform  distribution  with  depth,  he  was  able  to  match  a  number  of  detailed 
P-models;  however,  the  corresponding  predicted  S-modelr  bear  little  resemblance 
to  shear  models  determined  by  seismic  observations  (see  Figure  1).  Of  course, 
this  shear  profile  can  be  easily  altered  in  a  seemingly  infinite  number  of 
ways,  changing  the  Fe  content  as  a  function  of  depth  being  an  obvious  example. 
It  appears  that  much  of  this  ambiguity  can  be  eliminate.;  by  Inverting  detailed 
observed  profiles  of  P  and  S  together.  For  this  reason,  we  think  a  serious 
-eview  of  existing  shear  models  is  in  order. 

There  are  essentially  two  types  of  seismic  observations  used  in 
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shear  nodel  determinations,  namely,  eigenperiods  from  free  oscillations 
and  travel-times  of  body  phases.  Earth  structure  from  free  oscillations 
has  studied  extensively  (see  for  example  Jordan  and  Anderson.  1973, 

and  i..:iewonski  and  Gilbert.  19/2].  These  methods  are  quite  successful' 
in  determining  the  gross  structure  of  the  earth  but  probably  lack,  the 
resolving  power  necessary  for  detailed  upper  mantle  modeling.  Shear  models 
based  on  body  phases  are  relatively  rare,  some  major  contributors  are  Kovach 
and  Tohinson  [1969],  Robinson  and  Kovach  [1971],  Ibrahim  and  Nuttli  [1967] 
and  •Cuttli  [1972]. 

In  the  Kovach  and  Robinson  [1969]  study  shear  waveforms  are  not  used 
except  to  measure  the  apparent  velocity  of  the  starting  pulse  across  the  TFO 
LRSM  array.  We  repeated  many  of  their  measurements  and  found  good  agreement, 
but  it  should  be  realized  that  dt/dA  values  of  only  the  first  arrivals  does 
not  lead  to  unique  models  [for  example,  see  STAN2  of  Kovach  and  Robinson  [1969], 
or  Wiggins  [1970].  Attempts  at  measuring  the  dT/dA  values  for  second  arrivals 
were  not  particularly  successful  as  pointed  out  by  the  above  authors. 

On  the  other  hand,  Ibrahim  and  Nuttli  [1967]  and  Nuttli  [1972] 
used  particle  motion  diagrams  to  determine  the  onset  of  the  multi-arrivals 
produced  by  the  various  triplications.  The  method  Is  successful  in 
identifying  signals,  however,  it  appears,  as  demonstrated  later,  that  many 
such  arrivals  are  not  true-  jpper  mantle  shear  arrivals  but  are  PL-coupled 
shear  waves  and  other  P-SV  crustal  arrivals.  We  will  avoid  these  problems 
by  working  exclusively  with  the  SH  component  although  this  creates  considerable 
observational  difficulties. 


SYNTHESIZING  BODY  HAVES 


In  order  to  interpret  the  rather  complicated  seismograms  we  must 
separate  the  various  contributions  due  to  upper  mantle  structure  from 
other  effects  such  as  source  complexity.  To  facilitate  this  task  we 
calculate  synthetic  seismograms  assuming  various  types  of  sources.  The 
problem  can  be  somewhat  simplified  by  computing  the  SH  component  and 
performing  the  synthesizing  in  two  stages.  First,  we  compute  a  step 
function  response  of  a  layered  model  neglecting  the  radiation  field  and 
surface  effects  [see  Helmberger  [1973a]  and  Gilbert  and  Helmberger  [197]]. 
Secondly,  we  generate  an  effective  source  function  which  Includes  the 
surface  interaction  for  an  event  at  a  given  depth  and  vertical  radiation 
field  by  assuming  a  particular  ray  parameter  as  discussed  by  Helmberger  [1974]. 
This  procedure  allows  us  to  store  the  relatively  expensive  step  function 
response  and  test  various  source  configurations  by  convolution.  Note  that 
this  procedure  is  not  particularly  accurate  fov  cases  where  the  vertical 
radiation  pattern  is  varying  rapidly  at  ray  parameters  appropriate  for 
upper  mantle  triplications.  He  will  attempt  to  avoid  such  sources  in 
this  study  by  relying  mostly  on  strike-slip  events.  The  effective 
source  for  the  strike-slip  and  dip-slip  cases  as  a  function  of  epicentral 
depth  is  given  in  Figure  2.  The  synthetics  for  depth  ■  6  km,  SH6,  are 
appropriate  for  the  Borrego  Mountain  Event  [see  Helmberger .  1974]. 

Synthetic  responses  for  models  STAN3,  IN,  and  US26  are  given  in  Figures  3,  4, 
and  5  using  SH6  as  the  effective  source. 

The  general  features  of  these  profiles  are  controlled  by  the 
relative  Jumps  in  velocity  at  the  two  transition  zones.  In  Figure  4 
the  response  is  almost  completely  dominated  by  the  two  generalized  rays 
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vK’scribing  this  Interaction.  The  reflections  are  Indicated  by  the  peaks 
‘  r*.  •  refractions  by  the  preceding  ramps.  These  two  responses  cross  near 
i-v  .  Fn  a i-\N3,  the  transition  zones  are  much  more  subdued  and,  consequently, 
the  reflected  pulses  are  less  obvious.  The  responses  cross  near  20*  and 
produce  a  rather  strong  synthetic,  a  very  Important  feature  as  we  will  see 
shortly.  The  US26  model  is  a  reinterpretation  of  the  travel - 
time  data  of  Ibrahim  and  Nuttll  [1967J  to  produce  a  shear  model  with 
transition  zones  compatible  with  Johnson’s  [1967]  P  velocity  structure. 

The  low  velocity  zone  starting  at  500  km  Is  required  to  fit  the  triplication 
data.  The  synthetics  from  this  model  are  similar  to  those  from  model  IN. 

Models  STAN3  and  IN  although  derived  from  essentially  the  same  data  bank, 
differ  In  travel  times,  triplication  points  and  waveshape.  In  fact, 
their  predicted  waveshapes  are  remarkably  discordant. 

OBSERVATIONS 

The  observations  used  In  this  study  are  the  long-period  LRSM  and  WWSS 
responses  to  earthquakes  and  nuclear  events.  Plgures  6  and  7  show  the  locations 
of  the  various  stations  and  events.  Moec  of  the  data  Is  relevant  to 
'.Vs tern  United  States  although  we  Included  the  Arctic  event  to  sample 
some  of  the  expected  lateral  variation.  Depth  coverage  is  from  about 
150  to  1500  km.  The  region  sampled  is  predominantly  that  used 
in  the  P-wave  inversion  discussed  by  Helmberger  and  Wiggins  [1971],  and 
'■•'i  gains  and  Helmberger  [1973]. 

Samples  of  the  WWSS  observations  used  are  given  In  Figures  8,  9  and 
10.  Some  of  the  data  In  Figure  9  has  been  rotated  Into  SH  whereas  In 
Figures  8  and  10  the  observations  are  predominantly  SH  by  choice  of 
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station  component.  Complete  seismo?r<*ms  of  some  sample  LRSM  observations, 

.  the  SH  pulse  emerging  from  the  Love  wave  are  displayed  in 

Figures  11  and  12.  Before  giving  a  detailed  account  of  the 

data  we  will  first  discuss  the  arguments  against  using  SV  type  motions. 

Due  to  a  fortunate  set  of  circumstances  the  observations  from  the 
Borrego  Mountain  earthquake  as  recorded  in  southeastern  Camda  are 
naturally  polarized  into  SH(NS)  and  SV(EU)  (Figure  13).  As  expected  from 
knowledge  of  the  radiation  pattern  of  this  strike  slip  event,  we  find  the 
first  downgoing  peak  of  the  SV  component  to  be  roughly  half  the  amplitude 
of  SH  at  SCB  and  slightly  less  at  STJ.  There  is  good  agreement  between  the 
SH  pulses  and  the  synthetics  in  Figure  2,  SH6.  The  glich  occurring  in  the 
1st  peak  can  be  interpreted  as  the  i  flection,  sS,  from  the  surface 
[Helmberger,  1974].  Unlike  the  SH  pulses  which  do  not  appear  to  change 
shape  with  distances,  the  SV  waveform  is  elongated  and  shows  a  complex 
behavior  as  a  function  of  distance,  becoming  somewhat  simpler  at  A  ■  49°. 

V.’i  interpret  this  interesting  behavior  as  primarily  due  to  shear-coupled 
PL  waves  and  not  the  result  of  an  upper  mantle  triplication  as  previous 
workers  have  suggested. 

Extensive  literature  exists  on  the  theory  of  PL  waves  [see  for 
example  Phinney,  1961;  Su  and  Dorman,  1965].  Most  PL  waves  are  produced 
by  the  coupling  between  P  and  SV  waves  in  the  crustal  waveguide.  The 
type  discussed  have  followed  the  SV  wave,  displaying  prograde  particle 
motion,  and  are  sometimes  called  shear-coupled  PL  waves  a*ter  Chander  et^  al . 
[1968].  Particle  motion  diagrams  along  with  the  product  trace  (Z*EV)  used 
by  Ibrahim  and  Muttll  [1967]  and  Robinson  and  Kovach  [1971]  for  arrival 
identl f icatlona  are  displayed  in  Figure  13.  The  onset  of  the  PL  phase  can  be 
determined  by  the  start  of  the  prograda  motion  which  is  especially  easy  to 
Identify  on  the  SCB  record.  Arrows  mark  the  beginning  of  PL  on  the 
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various  traces.  Note  the  change  of  period  as  a  function  of  distance.  This 
type  of  dispersion  has  been  studied  by  Chander  et  al .  [1968]  and  can  be  used 

t 

to  establish  the  range  at  which  the  PL-coupling  occurs.  Uslpg  their  figure 
on  phase  velocity  versus  period  indicates  that  for  the  SCB  observation  the 
.  . '  in.  occurred  near  20°.  The  PL  onset  at  this  range  is  about  20  sec  behind 
Ir-,:  arrival  and  falls  on  the  travel  time  branch  associated  with  the 
bottom  transition  zone  as  defined  by  Ibrahim  and  Nuttll  (1967),  thus  the  confusion. 
In  fact,  a  number  of  workers  have  pointed  out  that  for  some  reason  the  arrival 
from  the  bottom  transition  is  most  easily  identified  on  the  vertical  or  radial 
component,  for  example  see  Figure  8  of  Kovach  and  hobinson  [1969]. 

There  are  other  difficulties  with  working  with  SV  waveshapes  such  as 
converted  P  waves  arriving  before  the  direct  SV  [see  Kanasewlch  et  al.,  1974]. 

This  feature  can  be  seen  in  Figure  13  where  the  S  onset  is  unclear.  There 
is  also  the  problem  of  waveform  distortion  caused  by  receiver  function  complexity 
near  203 .  That  is,  when  the  reciprocal  of  the  SV  ray  parameter  approaches 
the  surface  compressional  velocity  one  obtains  a  complex  interaction.  The 
result  is  a  waveshape  distortion  that  ia  a  function  of  receiver  geology  and, 
obviously,  difficult  to  assess.  For  the  above  reasons  we  will  rely  entirely 
on  SH  although  this  made  the  data  search  much  more  demanding. 

Table  1  lists  the  travel  'imes  and  events  used  in  this  study.  The 
times  have  been  corrected  for  source  depth  and  instrumental  delay.  That  is, 
we  found  tnat  the  WVJSS  instrument  has  about  a  two  second  lag  in  its  transient 
response  compared  to  three  seconds  for  the  LRSM  instrument.  These  values 
were  determined  from  Figure  2  so  they  also  contain  the  correction  for  an  assumed 
Q  and  source  description.  The  data  are  plotted  in  Figure  14.  We  included 
t  ‘  il  data  from  NTS  events  as  given  by  Nuttli  [1969]  and  found  good  agreement. 

The  data  from  the  Arctic  event  plot  near  the  Borrego  observat  ion.*,  it 
large  ranges  but  plot  progressively  early  moving  towards  smaller  ranges 
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. -aro  14).  This  feature  suggests  shallow  lateral  vari.  tions  as  might 


1  by  the  absence  of  a  low  velocity  zone.  This  conclusion  is 


instantiated  by  the  presence  of  a  very  strong  Sn  or  high  frequency  S 


:  r-’nr.es  less  than  20°  for  shield  observations.  The  waveforms  beyond 


iri  dominated  by  deep  structure  and  appear  similar  to  the  other 


'starvations. 


SHEAR  MODEL 


Our  search  technique  for  finding  acceptable  models  was  essentially 
trial  and  error  with  the  goals  of  fitting  the  travel  times  of  Figure  14 


ind  observed  high  quality  SH  waveforms.  We  were  guided  in  this  search  by 


'•imi lari  ties  in  previously  generated  synthetic  seismograms  and  observed 
waveforms.  For  instance,  the  synthetic  at  A  -  14°  for  model  STAN 3  (Figure  3) 


looks  similar  to  the  observations  at  A  -  15. 6*  displayed  in  Figure  12.  From 


this  comparison,  one  can  roughly  determine  the  relative  position  of  the 


triplications  by  examining  the  step  function  response  of  S'  1.  Proceeding 


in  this  fashion  we  obtained  the  models  given  in  Figure  15.  The  corresponding 


travel  time  triplications  for  the  preferred  model  SHR14  are  given  in 


Figure  14  for  comparison  with  the  first  arrivals.  The  synthetics  are 


displayed  in  Figures  16  and  17. 


Some  of  the  general  features  of  our  final  model  SHR14  are  quite 


''  r  10  those  of  STAN3.  The  two  major  triplications  associated  with 


the  400  and  650  km  transition  zones  still  exist  although  considerably 


subdued.  Reflections  from  these  two  zones  cross  at  20°.  The  small 


triplication  produced  by  the  velocity  Jump  near  500  km  in  conjunction 


with  the  constructive  interference  caused  by  the  above  cross-over 


produces  extremely  strong  shear  waves  near  20*.  Note  the  amplitudes 
near  20  are  nearly  an  order  of  magnitude  greater  than  at  30*. 


I 
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-0"'  the  synthctlcs  *re  dominated  by  the  first  arrival  with 
son.  interference  near  28”  and  with  little  aign  of  a  second  arrival 
beyond  30° . 

The  major  features  of  the  data  that  dictate  this  model  are  most 
easily  discussed  with  respect  to  Figure  14,  where  the  major  triplication 
points  and  branch  crossings  are  labelled.  By  fitting  the  observed 
waveforms  w.  verify  the  relative  position,  of  these  critical  points  and 

associated  model..  We  will  devote  the  remainder  of  this  section  to  such 
a  detailed  comparison. 

Comparing  the  .ynthetlca  generated  in  Figure  16  with  observations 
requires  a  great  deal  of  patch  work  due  to  the  lack  of  a  complete  profile 
of  da,,.  This  situation  requires  using  short  profiles  of  observations 
"ic  articular  source  functions;  one  such  comparison  is  given  in 
Figure  18.  By  fitting  these  observation,  we  hope  to  determine  the  position 
of  branch  DB  with  respect  to  branch  AB.  This  determination  is  especially 
important  since  the  apparent  velocities  as  measured  by  first  arrivals 
between  BE  are  considerably  higher  than  the  branch  DB  indicates.  First, 
we  note  that  the  two  LRSM  observations  in  Figure  11  are  very  much  alike, 
and  that  they  also  agree  with  the  synthetic,  of  Figure  16.  We  assume 
that  our  SH6  source  function  1.  appropriate  and  group  these  two  events 
together  at  shorter  ranges  a.  given  in  Figure  18.  The  dashed  line  in 
this  figure  is  the  same  as  branch  AB  in  be  travel  time  plot.  The  scatter 
in  time  picks  as  given  in  that  plot  can  be  easily  seen  by  examining  the 
observation,  in  Figure  18.  There  are  alao  con.ider.ble  differences  in 
waveshape,  although  .ome  „f  thI.  ,ffect  prob.bly  In8truMnul  „  Jn  £he 

case  of  HRAZ.  The.,  dl.tortlon.  or  vave.hap.  in.t.bllltie.  could  be 


- 
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C”US<!li  by  lateral  Structure  °r  perhaps  vertical  radiation  at  the 
-otto.  We  Win  he  primarily  concerned  with  the  first  50  sec.  of  waveform 
nee  the  latter  portion  of  the  recordings  are  Influenced  by  shallow 
lecture  fLove  waves,  see  Figure  12).  it  „  rather  difficult  to  establish 
.  ..  T.ooiness  of  fit  of  the  synthetics  to  the  observations  without  a  formal 
on,....  on  hut  we  will  character Ire  the  fit  by  the  relative  separation  between 

tlU'  ""i0r  i,05U1Ve  P“'“  SS  “  °f  distance.  A  traced  overlay  is  useful 

■or  this  purpose.  The  fit  Is  good  although  the  synthetic  first  arrival  at 

16°  appears  small.  Parr  nf  rm.  aa 

of  this  disagreement  is  caused  by  Inadequate  ray 

rno  e*in8  in  that  we  neglected  the  vertical 

rtlcal  radiation  pattern  which  would  enhance 

the  first  arrival  about  102. 

The  waveshapes  of  the  UVSS  observations  are  considerably  more 
—sting  than  the  LHSM  example,  of  Figure  18  although  they  are  rather 
sparse  (see  Figure.  9  s„d  10).  Most  of  these  event,  are  strike-slip 
determined  In  the  literature  and  can  be  compared  directly  with  the  synthetics 
of  Figure  16.  However,  their  source  time  histories  are  expected  to  be 
somewhat  different  due  to  variation.  In  source  depth  and  dislocation 
history.  The  match  of  the  synthetic  to  observed  shape,  of  Figure  10  Is 
particularly  convincing.  The  Interference  In  the  first  peak  of  the 
synthetic  near  17-  is  observed  a.  well  a.  the  apparent  pulse  broadening 
^Should  also  be  noted  thaf  the  SH  amplitudes  increase  rapidly  approach,,, 

0  “  Can  b*  See"  by  co"Parl"8  ^e  SH  pulse  directly  with  the  Love  wave 
It  Is  not  uncommon  to  observe  the  SH  weve  stronger  than  the 
trace  waves  near  20*.  This  strong  focu.in,  „f  energy  .How.  relatively 
small  earth, u.k.a,  magnitude.  1...  than  5,  to  be  examined  «  this  r.nge. 
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,0r  •'>:"apl°  ”*  the  "■  «•>•*«•“«  of  -gnitudc  5  at  21.4®  ,n  Flgure 
I  1 1 tud(  s  near  20  can  also  be  observed  In  Nuttll  (1969)  data 
<-•  His  Figure  6).  The  synthetics  fro.  .ode!  SHB14  dlspl.^  these  features 
V* TV  veil.  This  is  not  the  case  for  the  models  IN  and  US26.  STAN  3  predicts 

strong  amplitudes  near  20*  but  the  shape,  sre  somewhat  discordant, 
especially  at  17*. 

Another  comparison  of  LKSM  observations  with  synthetic.  Is  given 
in  Figure  19.  We  assumed  s  dlp-,l,p  .ource  DS14  (see  Figure  2)  since  this 
effective  source  looked  the  Boat  like  the  observed  record  st  21*  s„d  1. 
compatible  with  the  F  wave  rsdlstion.  We  explain  the  elongated  signal  st 
26»  as  a  combination  of  branches  BC  and  CD. 

Beyond  20®  It  become,  much  more  difficult  to  obt.ln  quality  dst. 

<iuc  to  the  shear  coupled  PL  w.ve  ..  discussed  earlier.  W,  have  to  rely 

to  resolve  SH  sni  SV  since  rotating  does  not 
t0°  CffeCt‘V'-  Tl"  portion  of  the  BLC  waveshape  (Figure  9, 

is  dotted  because  of  possible  cont.mln.tion  by  strong  SV  motion  (PL  problem) 

The  front  portion  of  the  recording  Indicting  interference  appears  accurate 
Since  this  shape  Is  reproduced  by  many  events  such  a.  the  observed  shape 
ut  AAM  from  the  Borrego  event.  The  Interference  at  this  range  Is  the 
'I.  mutation  of  the  presence  of  a  triplication.  This  feature  Is  used 
• •  ■  '  V  branch  HI.  Faced  with  the  lack  of  quality  waveshapes  at 

,h'M  raaR,’S  UO  a”°"Pt-d  t0  •«*“«  domain  of  possible  models  that 
■atis.y  the  travel  time  data,  the  28®  Interference,  snd  waveshape,  beyond 

'0‘  <l3Ck  °f  SeC°nd  ‘rrlval'>-  "»  “del  ED42  given  Pigure  „  „ 
attempt  at  Including  a  low  velocity  ton.  between  400  and  600  km.  The 
interference  at  28®  1.  relatively  es.y  to  produce!  however,  to  fit 
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’In  > •  i'  i  1  times  requires  that  the  600  kra  transition  be  very  large 

whi  h  produces  strong  second  arrivals  beyond  30°.  One  can  suppress  this 

feature  bv  making  the  first  arrival  stronger  (see  Figure  17).  Yhe  Increase 

in  velocity  occurring  at  850  km  produces  this  effect  from  28  to  32°  but  the 

problem  recurs  at  35°  (see  Figure  20).  Synthetics  for  some  other  models 

are  included  for  comparison.  We  don't  feel  that  the  data  supports  much 

of  a  second  arrival  at  these  ranges  (see  Figure  8) .  This  model 

also  predicts  complicated  waveforms  from  20  to  24°  as  does  the  models 

l*S26  and  IN.  The  first  arrival  is  essentially  diffracted  energy 

and  is  thus  small  compared  to  the  strong  arrival  returning  from  the 

600  km  transition. 

On  the  other  hand,  most  observations  in  the  range  20  to  24°  appear  quite 
simple,  for  example,  see  MNN  in  Figure  9.  This  observation  is  especially 
clear  since  it  occurs  on  a  P-SV  radiation  node.  There  is  littxe  evidence 
of  any  motion  on  all  three  components  other  than  SH  and  Love.  However,  a 
small  first  arrival  could  still  be  lost  in  the  background  noise.  To  resolve 
the  question  we  examined  relatively  strong  signals  as  displayed  in  Figure  21. 

The  OXF  observation  of  the  Borrego  event  which  is  well  located  has  a  back 
azimuth  of  274°.  The  trace  goes  off  scale  as  expected  but  there  does  not 
appear  to  be  a  significant  precursor.  The  same  is  true  for  the  5  July  event 
as  recorded  at  MNN.  This  event  is  also  on  a  P-SV  radiation  node.  A  detailed 
-  ire':,  or  this  small  arrival  proved  fruitless.  For  this  reason,  we  conclude 
that  the  strong  arrival  occurring  within  the  first  few  seconds  of  onset  at  these 
ranges  is  coming  from  the  600  ka  transition  and  that,  therefore,  the  velocity 
must  Increase  below  500  km.  Any  decrease  would  move  this  strong  arrival 
back  in  time  as  demonstrated  in  Figure  17. 
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DISCUSSION 


At  this  point  it  seems  useful  to  compare  the  travel  times  predicted 
by  the  proposed  models  and  discuss  the  various  disagreements  with  respect 
to  inadequate  inversion  techniques  versus  problems  of  lateral  variation. 

The  reduced  travel  times  relative  to  the  Jeff reys-Bullen  [ 19581  tables  are 
displayed  in  Figure  22  assuming  a  surface  focus.  The  difference  between  IN 
ii  4  in  travel  times  is  due  to  interpretation  of  seismograms  as 

t  earlier.  There  does  not  appear  to  be  any  significant  difference 

in  si!  travel  times  as  given  by  Nuttli  [1969]  even  though  we  used  different 
events .  We  think  the  reason  for  this  is  that  he  used  explosions  which  are  well 
located  and  we  used  mainly  well  located  and  timed  earthquakes  such  as 
Borrego  Mountain,  San  Fernando,  and  Truckee.  The  travel  times  predicted 
from  the  STAN 3  model  are  substantially  faster.  However,  the  change  of 
slope  was  the  actual  measurement  and  the  rapid  increases  at  20  and 
J4C  are  also  apparent  in  SHR14.  The  difference  between  17°  and  19°  is 
due  to  the  disagreement  in  branch  crossings  at  B  (see  Figure  14).  The 
above  differences  are  mostly  disagreements  in  data  handling  and  probably 
not  caused  by  lateral  variation.  On  the  other  hand,  the  travel  time 
difference  between  SHR14  and  the  Jef freys-Bullen  [1958]  tables  could 
easily  be  caused  by  regional  structure.  That  is,  the  travel  times  appropriate 
for  midwestern  United  States  may  be  anomalous  compared  to  Europe  and  other 
continental  shield  regions  as  suggested  by  the  results  from  the  Arctic 
event  (Figures  8  and  14).  There  is  also  the  problem  of  distinguishing 
between  Sn  and  S  when  processing  shield  observations  for  ranges  less  than 
20°.  This  may  be  the  reason  why  the  major  transition  zone  at  400  km  was 


not  discovered  earlier. 
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RESULTS 


14. 


We  have  attempted  to  determine  a  shear  model  that  fits  both 
travel  times  and  waveshapes,  essentially  by  trial  and  error.  Fittine  fine 
->(  rur'  i- •  hv  this  method  is  extremely  difficult  and  somewhat  inefficient.  '.'or 
tli-  mi  we  have  tried  to  produce  a  starting  model  for  formal  inversion 

i  iw  being  developed  by  exploring  the  possible  domains  of  solur ions 
.-iv., i  1 .1  r  lc.  We  think  SHR1A  is  a  promising  candidate.  A  slightlv  smr-'  ed 
v.  r .  ,1  long  with  a  P-model,  essentially  L1A  proposed  by  Helmberger  and 
-  i gg ins  [1971]  that  has  transition  zones  in  roughly  the  same  locations, 
t.s  given  in  Figure  23.  We  also  included  a  mod  1  proposed  by  Jordan  aid 
_rderson  [1967]  for  comparison.  The  model  parameters  are  given  in  Table  2. 

In  summary,  we  will  review  the  major  structural  features  of  our  model 

that  predict  realistic  seismograms  fitting  the  data.  However,  due  to  the 

data  distribution  we  have  determined  some  parts  of  the  model  considerably 

better  than  others.  We  will  attempt  to  convey  the  features  which  ve  think 

are  well  constrained  by  this  study.  First,  the  structure  near  the  surface 

is  taken  from  Stewart  and  Pakiser  [1962]  who  have  observed  both  Pn(8.2,  and 

S  (A.A5)  on  the  same  profile  of  the  Gnome  explosion.  These  velocities  are 
n 

in  agreement  with  measurements  made  in  southern  California  using  quarry  blasts 
presently  being  investigated.  The  structure  of  the  low  velocity  zone,  LVZ, 
taker,  from  Helmberger  [1973].  But  since  the  structure  of  the  LVZ  is  highly 
variable  in  western  United  States  [York  and  Helmberger.  1973],  this  model 
can  only  be  considered  as  representative.  Fortunately,  the  details  of  the  LVZ 
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'o  not  greatly  influence  our  upper  mantle  structure  except  by  slight  vertical 
shifting.  The  velocity  below  200  kra  is  fixed  by  the  (dT/dA)  values,  essentially 
the  slope  of  the  branch  AB  in  Figure  14.  This  branch  appears  to  start  at 
about  13J,  in  that  the  first  arrival  contains  higher  frequency  than  the  correspond¬ 
ing  Love  wave  with  these  high  phase  velocities  (see  Figure  9) .  At  larger 
ranges,  the  AB  branch  bends  downward  indicating  a  positive  velocity  gradient. 

This  feature  suggests  relatively  strong  first  arrivals  which  can  be  seen 
in  Figure  18.  The  overall  fit  of  the  synthetics  to  observed  waveforms  over 
.  «.•  range  14  to  19°  is  quite  good  which  leads  us  to  conclude  that  the  model 
fron  200  to  450  km  is  reasonably  accurate.  (We  will  postpone  questions  about 
hew  accurate  to  a  later  date.)  The  transition  at  500  km  is  based  mostly 
on  amplitude  and  waveshape  data.  Since  the  velocity  must  increase  substantially 
below  450  km  to  obey  the  first  arrival  times  (s  =e  Figure  14)  we  simply  added 
this  structure  to  produce  the  proper  waveshapes  and  produce  the  large  amplitudes 
near  20°.  At  larger  ranges  the  apparent  velocity  increases  rather  abruptly 
i.  e  Figure  14)  at  about  23°.  Similar  values  were  obtained  by  Kovach  and 
. ob inson  [1969]  as  is  evident  from  the  slope  of  the  reduced  travel  time  curve 
displayed  in  Figure  22.  The  magnitude  of  the  600  km  transition  zone  can 
be  determined  by  the  jump  in  velocity  whereas  its  smoothness  can  be  ascertained 
by  the  lack  of  prominent  second  arrivals  as  Indicated  in  Figure  8.  At  greater 
a  pths  the  velocity  approaches  Bl. 

In  conclusion,  a  preliminary  model  of  both  compressional  and 
shear  velocities  is  presented.  The  importance  of  the  seismic  parameter, 


a 
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16. 

determined  from  these  velocities,  In  assessing  phase  changes  versus 
compos  1 1 lonal  changes  is  obvious.  We  expect  the  model  to  be  elcered 
sotk wliat  with  the  addition  of  wide-band  array  measurements  and  better 
inversion  techniques.  Nevertheless,  a  mineralogic  interpretation  of 
these  combined  profiles  should  prove  interesting. 
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Table  1. 


Location 

I.at.  Long. 

Origin  Tine 
h  ra  a 

Station 

Distance 

A, deg 

Travel 

(min 

Time 
&  sec) 

Al.1’  '3  J 

4  ..i  N 

126.2  W 

09/27/9.3 

MNN 

24.05 

9’ 

*.5" 

MDS 

26.69 

10 

21" 

31  Mar  o4 

50.8  N 

130.2  W 

09/01/30 

GOL 

20.60 

8' 

37" 

MDS 

28.63 

10’ 

42" 

5  Jui\  64 

26.2  N 

110.2  W 

19/7/59 

BOZ 

19.25 

8' 

8.5 

COR 

21.05 

8' 

44" 

LON 

22.35 

9' 

11" 

MNN 

23.04 

9' 

19" 

ATL 

23.40 

9 

28" 

l"l  July  64 

44.7  N 

129.9  W 

06/47/54 

MNN 

24.25 

9’ 

40" 

j  o  •  >  ■  *4 

43.5  N 

127.5  W 

13/59/37 

KNUT 

12.94 

5' 

34" 

SGAZ 

13.49 

5' 

46" 

JRAZ 

14.81 

6' 

22" 

LGAZ 

15.37 

6’ 

35" 

HRAZ 

15.69 

6’ 

45" 

GEAZ 

16.42 

7’ 

1" 

RTNM 

18.92 

7’ 

59" 

1  Our  64 

45.7  N 

122.8  W 

12/31/24.6 

JRAZ 

14.44 

6’ 

11" 

LGAZ 

14.99 

6’ 

25" 

GEAZ 

16.04 

6' 

51" 

RTNM 

18.49 

V 

50" 

LUB 


21.87 


8'  59" 
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Table  1  (continued) 


'  lie 

Location 

Lat.  Long. 

Origin  Tine 
b  ra  a 

■ 

31.9  N  117.1  W 

20/54/35.3 

Sct't  66 

iO.  4  N 

120.1  W 

16/41/1.7 

Nov  67 

80.2  N 

1.  W 

13/42/02 

Station 

Distance 

4, deg 

Travel 
(min  & 

Time 

sec) 

MNN 

22.70 

9' 

12" 

FLO 

22.81 

9’ 

15" 

MDS 

24.38 

9’ 

40" 

FSJ 

15.31 

6' 

40" 

JCT 

18.80 

7’ 

59" 

FFC 

19.60 

8' 

4" 

MDS 

23.04 

9' 

19" 

UME 

17.48 

7’ 

17" 

KON 

20.90 

8' 

33" 

NUR 

21.16 

8’ 

39" 

FBC 

24.29 

9' 

36" 

CMC 

27.75 

10' 

31" 

BLC 

28.33 

10' 

43" 

YKC 

32.78 

11' 

53" 

LOR 

33.10 

11' 

58" 

FCC 

33.27 

11* 

57" 

STJ 

37.17 

13’ 

8" 

PTO 

39.29 

13’ 

2S" 

TOL 

40.46 

13' 

50" 

SES 

43.99 

14' 

38" 

PHC 

45.93 

15’ 

12" 

PNT 


46.27 


15’ 


16" 
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Table  1  (continued) 


No.  Date 


1  Dec  67 


Location 
Lat.  Long. 

44.2  N  128.8  W 


•i  Apr  68 


33.1  N  116.1  W 


Origin  Tine 

Station 

distance 

Travel 

Time 

has 

A, deg 

(min  & 

sec) 

06/26/16 

DUG 

12.52 

5’ 

24" 

SES 

13.56 

5’ 

48” 

EDM 

13.61 

5' 

50" 

GOL 

17.98 

7' 

36" 

TUC 

18.43 

7' 

50" 

ALQ 

19.5 

8' 

13" 

YKC 

20.12 

8' 

18" 

COL 

23.30 

9 ' 

22" 

LUB 

23.37 

9' 

27" 

CMC 

24.77 

9' 

45" 

02/28/59 

LON 

14.30 

6’ 

9" 

OXF 

22.21 

9' 

5" 

ATL 

26.51 

10’ 

13" 

AAM 

27.13 

10' 

22" 

BLA 

29.35 

10' 

59" 

SCB 

30.55 

11' 

21" 

SCP 

31.31 

11' 

36" 

OTT 

33.26 

12’ 

1" 

BLC 

33.62 

12' 

12.5 

OGD 

33.78 

12' 

12.5 

WES 

36.25 

12’ 

51.5 

SFA 

36.83 

12' 

56" 

SCH 

40.55 

13' 

50" 

1 


HAL 


41.81 


14’  12" 
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Table  1  (continued) 


Date 

Location 

Lit.  Long. 

OriRin  Tine 
h  m  a 

Station 

Distance 

Travel 
(min  & 

Time 

sec) 

:■  :,ov  7 n 

3.78  N  127  45  W 

3/11/43 

GOL 

16.96 

7' 

17V 

ALQ 

18.44 

7 ' 

49H 

JCT 

25.59 

10' 

2" 

FLO 

28.14 

10’ 

40” 

OXF 

30.69 

11’ 

24" 

9  Feb  71 

34.41  N  118.40  W 

17/00/42 

DAL 

18.06 

7’ 

38" 

FLO 

22.88 

9* 

15" 

SHA 

23.87 

9’ 

34.5 

.  .■  n 

49.6  N  129.5  W 

5/50/5.8 

LHC 

26.16 

10’ 

5" 
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Table  II 


Depth 

V 

P 

V 

8 

* 

0 

25 

6.80 

3.69 

28.09 

.291 

50 

8.20 

4.45 

40.84 

.291 

75 

7.81 

4.30 

36.34 

.282 

100 

/.92 

4.43 

36.56 

.272 

125 

8.02 

4.44 

38.04 

.279 

150 

8.14 

4.46 

39.74 

.285 

175 

8.24 

4.48 

41.14 

.290 

200 

8.32 

4.50 

42.22 

.293 

225 

8.36 

4.52 

42.65 

.293 

250 

8.41 

4.54 

43.25 

.294 

275 

8.45 

4.56 

43.68 

.294 

300 

8.49 

4.58 

44.11 

.294 

325 

8.53 

4.62 

44.30 

.292 

350 

8.57 

4.65 

44.62 

.291 

375 

8.63 

4.72 

44.77 

.286 

400 

9.10 

4.98 

49.74 

.286 

425 

9.53 

5.12 

55.87 

.297 

450 

9.57 

5.17 

55.95 

.293 

475 

9.62 

5.18 

56.77 

.295 

500 

9.71 

5.23 

57.81 

.295 

525 

9.92 

5.39 

59.67 

.290 

550 

10.07 

5.43 

62.09 

.295 

575 

10.13 

5.43 

63.31 

.293 

600 


10.19 


5.48 


63.80 


296 
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Table  II  (continued) 


Depth 

V 

P 

V 

s 

4 

0 

625 

10.29 

5.59 

64.22 

.290 

650 

10.41 

5.76 

64.13 

.279 

675 

10.74 

5.94 

68.30 

.279 

700 

10.93 

6.07 

70.34 

.277 

725 

10.98 

6.09 

71.11 

.277 

750 

11.02 

6.12 

71.50 

.277 

775 

11.06 

6.15 

71.90 

.276 

800 

11.11 

6.19 

72.35 

.274 

825 

11.15 

6.23 

72.57 

.273 

850 

11.19 

6.26 

72.97 

.272 

875 

11.23 

6.29 

73.36 

.271 

900 

11.27 

6.32 

73.76 

.270 

925 

11.31 

6.34 

74.32 

.270 

950 

11.35 

6.36 

74.89 

.271 

975 

11.39 

6.38 

75.46 

.271 

1000 

11.44 

6.40 

76.26 

.272 

1025 

11.49 

6.42 

77.07 

.273 

1060 

11.53 

6.43 

77.82 

.274 
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FIGURE  CAPTIONS 

Figure  1.  Comparison  of  models  proposed  from  shear  wave  observations 

with  a  predicted  shear  profile  based  on  theoretical  mineralogic 
considerations. 

Figure  2.  Effective  source  function  containing  surface  interaction  as  a 

function  of  epicentral  depth;  contains  time  history,  instrument 
and  absorption. 

rieure  3.  Synthetic  profile  based  on  STAN3  assuming  SH6  as  the  effective 
source  function. 

Figure  4.  Synthetic  profile  based  on  IN  assuming  SH6  as  the  effective 
source  function. 

Figure  5.  Synthetic  profile  based  on  US26  assuming  SH6  as  the  effective 
source  function. 

Figure  6.  Map  locating  stations  and  events.  Events  are  identified  in 
Table  1  according  to  date  (see  Table  1). 

Figure  7.  Map  locating  an  event  in  the  Greenland  Sea  relative  to  stations 
using  an  Arctic  projection. 

Figure  8.  Profile  of  observations  from  events  23  Nov  67  (Arctic)  and 
9  Apr  68  (Borrego) . 

Figure  9.  Sample  recording  of  WWSS  observations  of  simple  events.  Note 
that  FBC  and  BLC  are  of  the  Arctic  event. 

Figure  10.  Examples  of  large  amplitude  arrivals  near  20*. 

Figure  11.  LRSM  observations  of  two  simple  events  off  Oregon  as  recorded  by 


the  same  station. 
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[  •  !Tf  12. 

!  j  gure  13. 

Figure  14. 
Figure  15. 
Figure  16. 

Figure  17. 

Figure  18. 

Figure  19. 

Figure  20. 

Figure  21. 

Figure  22. 

■'igu.-c  .'3. 


Example  of  interference  caused  by  upper  mantle  triplication. 
Observations  are  of  one  of  the  events  displayed  in  Figure  11. 
Observations  of  the  Borrego  event  as  recorded  on  three  components. 
Particle  motions  as  well  as  a  component  product  Z*EW  is 
Included  for  studying  the  PL  problem. 

Reduced  travel  time  plot  based  on  SH  data. 

Proposed  upper  mantle  models. 

Synthetic  profile  based  on  SHR14  assuming  SH6  as  the  effective 
source  function. 

Synthetic  profile  baaed  on  ED42  assuming  SH6  as  the  effective 
source  function. 

Comparison  of  synthetics  with  observations  assuming  SH6  as  the 
effective  source  function. 

Comparison  of  synthetics  with  4  June  64  observations  assuming 
a  dip-slip  source  (DS14) . 

Comparison  of  proposed  models  at  35*  displaying  the  magnitude  of  the 
second  arrival. 

Example  seismograms  showing  the  strong  onset  of  arrivals.  Arrows 
mark  the  arrivals  as  determined  by  the  authors. 

Comparison  of  residuals  of  travel  times  as  referenced  to  the 
Jef frey-Bullen  tables. 

Comparison  of  SHR14  with  a  model  determined  from  free-osc i llat ions 
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